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MINIMIZATION  OF  THE  TRUNCATION  ERROR  BY  GRID  ADAPTATION 


NAIL  K.  YAMALEEV* 

Abstract.  A  new  grid  adaptation  strategy,  which  minimizes  the  truncation  error  of  a  pth-order  finite 
difference  approximation,  is  proposed.  The  main  idea  of  the  method  is  based  on  the  observation  that  the 
global  truncation  error  associated  with  discretization  on  nonuniform  meshes  can  be  minimized  if  the  interior 
grid  points  are  redistributed  in  an  optimal  sequence.  The  method  does  not  explicitly  require  the  truncation 
error  estimate  and  at  the  same  time,  it  allows  one  to  increase  the  design  order  of  approximation  by  one 
globally,  so  that  the  same  finite  difference  operator  reveals  superconvergence  properties  on  the  optimal  grid. 
Another  very  important  characteristic  of  the  method  is  that  if  the  differential  operator  and  the  metric 
coefficients  are  evaluated  identically  by  some  hybrid  approximation  the  single  optimal  grid  generator  can 
be  employed  in  the  entire  computational  domain  independently  of  points  where  the  hybrid  discretization 
switches  from  one  approximation  to  another.  Generalization  of  the  present  method  to  multiple  dimensions 
is  presented.  Numerical  calculations  of  several  one-dimensional  and  one  two-dimensional  test  examples 
demonstrate  the  performance  of  the  method  and  corroborate  the  theoretical  results. 

Key  words,  truncation  error,  grid  adaptation  criterion,  finite  difference  approximation,  error  equidis- 
tribution 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  Grid  adaptation  has  now  become  widespread  for  solving  multi-dimensional  partial 
differential  equations  in  arbitrary-shaped  domains.  One  of  the  most  important  problems  associated  with  the 
adaptive  grid  generation  is  an  essential  effect  of  the  grid  point  distribution  on  error  in  the  numerical  solution. 
Until  the  present  time  little  attention  has  been  paid  to  the  fact  that  the  concentration  of  grid  points  in  regions 
which  most  influence  the  accuracy  of  the  numerical  solution  may  at  the  same  time  introduce  additional  error 
due  to  the  grid  non-uniformity  [1]— [3] . 

There  are  two  basic  strategies  of  the  grid  adaptation,  namely,  grid  refinement  and  grid  redistribution 
techniques.  In  the  first  approach  grid  nodes  are  added  to  locally  enrich  the  grid  to  achieve  higher  accuracy. 
In  the  second  approach  the  number  of  grid  nodes  is  fixed  and  the  idea  is  to  adjust  the  position  of  grid  points 
to  improve  the  numerical  solution  accuracy.  In  spite  of  significant  distinctions,  for  both  methods  reliable 
and  efficient  grid  adaptation  criteria  are  needed. 

A  number  of  grid  adaptation  criteria  based  on  the  equidistribution  principle  have  been  developed.  As 
has  been  shown  in  [4],  the  grid  point  distribution  is  asymptotically  optimal  if  some  error  measure  is  equally 
distributed  over  the  field.  One  of  the  widely-used  approaches  is  to  redistribute  grid  points  in  accordance 
with  the  arc  length  and  the  local  curvature  of  the  solution  curve  [5],  [6].  This  kind  of  clustering  is  intended 
to  reduce  the  error  in  the  vicinity  of  strong  gradients  and  local  extrema  of  the  numerical  solution,  but  it 
does  not  necessarily  guarantee  improvement  in  the  accuracy  where  the  solution  is  smooth. 

Another  class  of  methods  is  based  on  equidistribution  or  minimization  of  the  local  truncation  error  or 
its  estimate  [7]— [10].  In  [7]  the  error  estimate  obtained  by  using  a  finite  difference  approximation  of  the 
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leading  truncation  error  term  is  equidistributed  by  the  grid  point  redistribution.  Klopfer  and  McRae  [8] 
solve  a  one-dimensional  shock-tube  problem  using  the  explicit  predictor-corrector  scheme  of  MacCormack 
on  a  grid  dynamically  adapted  to  the  solution.  The  error  estimate  is  the  leading  truncation  error  term 
of  the  differential  equations  transformed  to  the  computational  coordinates.  The  metric  coefficient  is  taken 
as  a  linear  function  of  the  smoothed  error  measure.  For  solving  a  second-order  two-point  boundary  value 
problem  with  a  centered  second-order  finite  difference  scheme  Denny  and  Landis  [9]  suggest  to  determine 
the  optimal  coordinate  mapping  so  that  the  entire  truncation  error  vanishes  at  all  grid  points.  However, 
this  grid  generator  concentrates  grid  nodes  where  the  solution  is  smooth  rather  than  near  steep  gradients. 
Thus,  the  error  reduction  occurs  in  regions  which  do  not  practically  affect  the  numerical  solution  accuracy. 
An  alternative  technique  is  employed  in  [10]  where  the  optimal  coordinate  transformation  is  constructed  as 
the  solution  of  a  constrained  parameter  optimization  problem  minimizing  a  measure  of  the  truncation  error. 
The  error  measure  used  is  a  finite  difference  evaluation  of  the  third  derivative  of  the  numerical  solution 
calculated  in  the  computational  space.  The  main  drawback  of  all  the  methods  mentioned  above  is  the  fact 
that  the  error  estimates  do  not  properly  take  into  account  that  part  of  the  truncation  error  which  is  caused 
by  the  nonuniform  grid  spacing.  Furthermore,  it  is  not  clear  how  to  extend  these  methods  to  more  general 
equations  and  discretizations  as  well  as  to  multiple  dimensions. 

A  grid  adaptation  procedure  equidistributing  an  error  estimate  of  the  numerical  solution  has  successfully 
been  used  in  [11]  to  reduce  simulation  error  in  such  integral  quantities  as  the  lift  or  drag.  This  error 
estimate  is  directly  related  to  the  local  residual  errors  of  the  primal  and  adjoint  solutions  of  the  Euler 
equations.  As  it  follows  from  the  numerical  results  presented  in  [11],  the  order  of  accuracy  of  the  integral 
outputs  increases  by  one  if  the  proposed  adaptation  strategy  is  employed.  Although,  this  approach  provides 
significant  improvement  in  the  accuracy  of  the  functional,  the  error  estimation  procedure  is  quite  expensive 
in  terms  of  computational  time  since  except  for  the  solution  of  the  primal  problem  it  is  needed  to  solve  the 
adjoint  Euler  equations  that  doubles  the  computational  efforts. 

The  formulation  of  an  adaptive  mesh  redistribution  algorithm  for  boundary  value  problems  in  one 
dimension  has  been  presented  in  [12].  The  analysis  uses  the  error  minimization  to  produce  an  optimal 
piecewise-polynomial  interpolant  in  a  given  norm  that  leads  to  the  development  of  a  family  of  grid  adaptation 
criteria.  Despite  the  fact  that  the  present  approach  works  well  in  one  dimension  this  error  equidistribution 
analysis  can  not  be  directly  extended  to  multiple  dimensions  [13]. 

In  [14]  and  [15]  the  finite  element  residual  is  applied  to  provide  a  criterion  for  determining  where  a 
finite  element  mesh  requires  refinement.  As  has  been  noted  in  [16]  for  hyperbolic  problems  with  non-smooth 
solutions  the  finite  element  residual  may  be  an  ineffective  error  estimator  since  for  such  problems  the  residual 
measured  in  the  L2  norm  diverges  whereas  the  numerical  solution  converges  in  this  norm.  The  problem  might 
be  overcome  if  the  divergence  of  the  residual  is  localized  to  the  area  of  non-smoothness  and  the  residual 
would  then  be  used  as  a  local  error  indicator.  However,  the  localization  of  discontinuities  becomes  a  very 
complicated  problem  in  multiple  dimensions. 

It  can  be  shown  that  the  truncation  error  of  any  differential  operator  obtained  on  a  nonuniform  grid 
consists  of  two  different  parts.  The  first  one,  which  always  exists  on  a  uniform  mesh,  is  due  to  the  ap¬ 
proximation  of  the  differential  operator  itself.  The  second  one  is  caused  by  the  contribution  to  the  error 
from  the  nonuniform  grid  spacing.  As  the  grid  is  locally  refined  or  redistributed  the  first  part  of  the  error 
decreases  while  the  second  part  may  considerably  increase  because  of  the  grid  non-uniformity.  All  of  the 
equidistribution  methods  mentioned  above  redistribute  grid  points  in  accordance  with  one  or  another  error 
estimate  obtained  on  a  non- adapted  grid,  but  in  doing  so  the  grid  adaptation  itself  introduces  additional 
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error  which  changes  the  error  distribution.  Therefore,  to  account  for  this  change  in  the  error  distribution 
the  grid  adaptation  procedure  based  on  the  error  equidistribution  strategy  should  be  repeated  iteratively 
until  the  error  estimate  norm  is  equally  distributed  over  the  field.  Note  that  for  moving  meshes  dynamically 
adapted  to  the  solution  the  iterative  procedure  should  be  done  at  each  time  step  to  get  the  optimal  mesh 
characterized  by  having  the  error  equidistributed  throughout  the  domain. 

The  main  objective  of  this  paper  is  to  construct  an  optimal  coordinate  transformation  so  that  the 
leading  truncation  error  term  of  an  arbitrary  pth-order  finite  difference  approximation  is  minimized  that 
provides  superconvergent  results  on  the  optimal  grid.  In  contrast  to  the  error  equidistribution  principle, 
for  the  present  technique  a  posteriori  error  estimate  is  not  explicitly  required.  Furthermore,  the  new  grid 
adaptation  criterion  allows  one  to  minimize  the  error  due  to  the  differential  operator  itself  and  the  error  owing 
to  the  evaluation  of  the  metric  coefficients  simultaneously.  Another  very  attractive  feature  of  the  present 
approach  is  its  applicability  to  hybrid  approximations  which  depend  on  some  basic  properties  of  the  solution 
such  as  a  flow  direction,  sonic  line  and  others.  If  the  metric  coefficients  are  evaluated  by  the  same  hybrid 
discretization  used  for  the  differential  operator,  the  new  grid  adaptation  criterion  remains  valid  in  the  whole 
computational  domain  regardless  to  points  where  the  hybrid  scheme  switches  from  one  approximation  to 
another.  Extension  of  the  new  adaptation  criterion  to  multiple  dimensions  is  presented.  Numerical  examples 
considered  illustrate  the  ability  of  the  method  and  corroborate  the  theoretical  analysis. 

2.  Grid  Adaptation  in  One  Dimension.  We  consider  the  truncation  error  of  the  first  derivative 
approximated  on  a  ID  nonuniform  grid.  Let  x  and  f  denote  the  physical  and  computational  coordinates, 
respectively.  Without  loss  of  generality  it  is  assumed  that  a  <  x  <b  and  0  <  £  <  1.  A  one  to  one  coordinate 
transformation  between  the  physical  and  the  computational  domains  is  given  by 


(2.1) 

II 

g 

where 

(2.2) 

x(0)  =  a 
a;(l)  =  b. 

It  is  assumed  that  the  above  mapping  is  not  singular  so  that  the  Jacobian  of  the  transformation  is  a  strictly 
positive  function,  i.e. 

(2.3)  xe>0,  V£  G  [0, 1]. 

The  nonuniform  grid  in  the  physical  space  is  obtained  as  images  of  nodes  of  a  uniform  mesh  in  the  compu¬ 
tational  domain 

(2.4)  =  2;  (£z) ,  ,  i  0,1,.../. 

Taking  into  account  the  coordinate  transformation  Eq.(2.1)  the  first  derivative  of  a  function  f(x)  with  respect 
to  x  can  be  written  as  follows 

(2-5)  f.  = 

To  construct  a  pth-order  approximation  of  fx  in  the  physical  domain  we  approximate  fa  and  x(  by  some 
pth-order  finite  difference  expressions  in  the  computational  domain 

i+h 

E  «J  fi 

(2-6)  M/.>  =  -£r - ■ 

Pmxm 

m~i — tni 
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where  xm  =  x(£m),  fi  -  f(£i)i  Lh  is  a  finite  difference  operator;  the  indexes  h ,  l2  and  m1,m2  as  well  as 
the  coefficients  an  and  (3m  depend  on  particular  approximations  used  for  evaluating  and  X(,  respectively. 
Henceforth,  we  shall  assume  that  the  functions  /(£)  and  x(£)  are  smooth  enough  so  that  all  derivatives 
needed  for  the  derivation  are  continuous  functions  on  £  G  [0, 1].  Expanding  the  nominator  and  denominator 
of  Eq.(2.6)  in  a  Taylor  series  with  respect  to  &  and  omitting  the  index  i  on  the  right  hand  side  yield 


E  on  ft  =  /«  +  C'/4P+1)  Ae  +  0(A£P+1) 


(2.7) 

l=i-h 

i+m-z 

+  c;4,,+1)A?  +  o<Ae’+1) 

m=i-mi 

where 

(P+i)  _  dp+lx 
t  <9£p+i  ’ 

,<«>  _«■*'/  i£=i 

Cf  and  Cl  are  constants  dependent  on  at  and  /3m,  respectively.  Substituting  Eq.(2.7)  into  Eq.(2.6)  and 
taking  into  account  that  x £  >  0,  V£  G  [0, 1]  one  can  write 


(2.8) 


LhUx)  = 


h  +  C/A^/<p+1) 

st(i+q?^Hp+1)) 


+  0(A^+1). 


Assuming  that  A£  is  chosen  to  be  sufficiently  small  so  that  A^\x(p+1) /xK\  «  1,  Eq.(2.7)  can  be  linearized 
as  follows 


(2.9)  Lh{fx)  =  ^  (/«  +  C'Ae4P+1))  (l  -  Cl  ^^P+1))  +  0( Ae+1). 


Note  that  the  error  introduced  by  the  linearization  is  of  the  order  of  0(A£2p+2).  Neglecting  higher  order 
terms  in  Eq.(2.9)  we  have 

Ap+i)  x^p+1^ 

(2.10)  Tp(x)  =  Lh{fx)  ~fx=  ClAe-1 - Cl Ae-^fi- 

xi  x( 

The  right  hand  side  of  Eq.(2.10)  is  the  leading  truncation  error  term.  Thus,  if  the  metric  coefficient  x( 
is  evaluated  numerically  as  in  Eq.(2.6)  the  asymptotic  truncation  error  of  any  pth-order  finite  difference 
approximation  consists  of  two  different  parts,  one  of  which  is  due  to  the  evaluation  of  fa  and  the  second  one 
is  caused  by  the  discretization  of  the  metric  coefficient  x^.  It  should  be  emphasized  that  any  grid  adaptation 
based  on  minimization  or  equidistribution  of  the  first  part  of  the  truncation  error  alone  is  not  sufficient  since 
the  second  part  of  the  truncation  error  may  drastically  increase  in  regions  where  a:(£)  rapidly  changes.  In 
other  words,  any  inconsistent  grid  adaptation  transfers  the  error  from  the  first  term  of  the  truncation  error 
to  the  second  one  and  vice  versa.  To  minimize  both  parts  of  the  truncation  error  simultaneously  we  impose 
the  following  restriction  on  the  coordinate  mapping  x(f),  V£  €  [0, 1] 

(2.11)  | Clf<p+1)x(  -  Clx{p+1)fz\  <  0( A04 

If  Eq.(2.11)  holds  the  asymptotic  order  of  approximation  of  Eq.(2.6)-(2.7)  on  the  optimal  grid  generated  by 
the  mapping  a;(£)  is  p  +  1  in  the  entire  computational  domain.  Replacing  the  inequality  sign  by  the  equality 
one  in  Eq.(2.11)  the  grid  adaptation  criterion  can  be  expressed  as 

(2.12)  C//'p+S  -  Clxf+1)f,  =  0(A()xl 
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Recall  that  the  coefficients  C/  and  C*  depend  on  the  particular  approximations  used  and  do  not  depend  on 
/(£)  and  x(£).  One  of  the  most  important  classes  of  approximation  is  a  consistent  approximation  when  the 
same  difference  operator  is  employed  to  evaluate  the  derivatives  /$  and  x$.  In  this  case  the  coefficients  Cj, 
and  Cp  are  identical  and  Eq.(2.12)  is  simplified  to 

(2.13)  flp+1)  -  fxx<£+1)  =  OiAQxt 
or  setting  the  right  hand  side  equal  to  zero  yields 

(2.14)  4P+S  -  =  °- 

There  are  several  advantages  of  such  a  simplification.  First  of  all,  the  use  of  the  same  difference  approx¬ 
imation  for  both  and  eliminates  the  fx  term  from  the  truncation  error  which  is  the  most  troublesome 
part  of  the  error  being  dependent  on  the  first  derivative  which  is  evaluated.  Actually,  let  us  represent  and 
in  terms  of  the  x  derivatives 


=  (w+x^)1 


f  =  x<f+i)fx+(p+ 1)4 


+  . . .  +  X‘ 


p+l  A  p+l) 


Note  that  the  binomial  theorem  can  not  be  used  to  expand  the  power  of  the  derivative  operator  in  the  above 
formula  since  d/d£  and  x^d/dx  do  not  commute.  Substituting  the  above  expressions  into  Eq.(2.10)  the 
leading  term.of  the  truncation  error  Tp(x)  can  be  Written  as  follows 


(2.16) 


Tp(x)  =  (l C}p  -  C*v) +  Cl  [(p+  \)xfh 


From  Eq.(2.16)  it  is  clear  that  if  Cl  ±  CTp  then  the  truncation  error  depends  on  the  first  derivative  fx  being 
approximated.  That  is  why  it  is  very  important  to  evaluate  the  metric  coefficient  by  the  same  difference 
approximation  used  for  f(.  It  should  be  noted,  that  if  xi  is  approximated  by  the  exact  analytical  expression 
or  any  finite  difference  formula  different  from  that  which  is  employed  to  calculate  it  gives  rise  to  the  fx 
term  in  the  truncation  error. 

Another  advantage  of  the  consistent  approximation  of  and  X(  is  that  the  single  optimal  grid  in  the 
sense  of  Eq.(2.14)  can  be  generated  for  hybrid  discretization,  when  the  coefficient  C/  may  implicitly  depend 
on  the  function  /(£)•  The  identical  numerical  approximation  of  and  /$  removes  the  dependence  of  the 
optimal  mapping  on  points  in  the  physical  domain  where  the  hybrid  scheme  switches  from  one  approximation 
to  another.  If  this  is  the  case  the  optimal  grid  point  distribution  depends  only  on  the  order  of  approximation 
and  is  completely  independent  of  the  particular  finite  difference  formula  used. 

As  has  already  been  mentioned,  Eq.(2.14)  is  a  grid  adaptation  criterion,  but  at  the  same  time  this 
equation  can  be  treated  as  a  grid  generation  equation.  To  provide  the  existence  of  the  solution  of  Eq.(2.14) 
it  is  assumed  that  f(  >  e  >  0,  V£  €  [0,1],  and  /(£)  €  Cp+1[0, 1],  It  can  easily  be  seen  that  x(£)  = 
Ci/(0  +  c2  is  the  solution  of  Eq.(2.14),  but  this  trivial  solution  is  not  appropriate  since  it  means  that  f(x) 
is  a  linear  function  of  x  in  the  physical  space.  Another  problem  associated  with  the  solution  of  Eq.(2.14) 
is  boundary  conditions.  Theoretically,  to  find  the  unique  solution  of  Eq.(2.14)  p+l  boundary  conditions 
should  be  imposed  while  only  two  boundary  conditions  Eq.(2.2)  are  available.  In  spite  on  the  abovementioned 
difficulties  the  optimal  grid  generation  problem  Eqs.(2.14),(2.2)  can  be  solved  analytically  for  very  important 
cases  p  =  1, 2  and  the  approximate  analytical  solutions  can  be  obtained  for  higher  order  discretizations  p  >  3. 
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2.1.  First-Order  Approximation,  p  =  1.  For  a  first-order  accurate  approximation  p  is  equal  to  one 
in  Eq.(2.13)  which  takes  the  form 

(2.17)  =  °(A6xf  • 


Using  the  following  expression  for  the  second  derivative 


r  „  fux£ 

JX*  -  -3 


Eq.(2.17)  written  in  the  physical  space  is  reduced  to 


(2.18) 


0(A£)' 


Integrating  Eq.(2.18)  and  taking  into  account  the  boundary  conditions  £(a)  =  0,  £(&)  -  1  yield 


J'  fxxd% 

(2.19)  «*)  -  ~b - • 

J'  fxx^X 
a 

However,  to  satisfy  Eq.(2.18)  the  following  restriction  should  be  imposed  on  fxx 

b 

(2.20)  J  fxxdx  =  0( AO- 

a 

Since  >  0  from  Eq.(2.18)  it  follows  that  fxx  >  0.  Consequently,  Eq.(2.20)  means  that  the  second 
derivative  fxx  has  to  be  of  the  order  of  0( A£)  for  all  x  €  [o,  i>].  In  other  words,  if  f(x)  is  an  essentially 
nonlinear  function,  so  that  Eq.(2.20)  is  not  satisfied,  it  is  impossible  to  increase  the  global  order  of  accuracy 
of  fx  by  the  grid  point  redistribution. 

2.2.  Second-Order  Approximation,  p  =  2.  If  both  /5  and  x(  are  evaluated  identically  by  a  second- 
order  accurate  formula  the  grid  adaptation  equation  Eq.(2.13)  written  for  p  =  2  becomes 


(2.21)  f^(x(  -  =  °(A0xl 

Let  us  transform  the  derivatives  in  Eq.(2.21)  from  the  computational  space  to  the  physical  space 

(2.22)  ^  —  fxX$ 

=  fxxxx |  +  &fxxx(xZ£  *b  fxxt((- 

Substituting  Eq.(2.22)  into  Eq.(2.21)  we  have 

(2.23)  fxxx  x\  +  3  f„xH  =  O(A0xl 
Using  the  following  expressions  for  the  metric  coefficient  and  its  derivative 

xt  =  i 
=  ~Ht 

and  assuming  that  fxx  /  0,  Vrr  £  [a,  6]  Eq.(2.23)  can  be  rewritten  as 


(2.24) 


Jxx  see  Jxi 
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Since  a  decrease  in  the  last  term  in  the  above  equation  increases  the  approximation  accuracy  we  neglect  the 
0(A£)  term  and  integrate  the  left  and  right  hand  sides  of  Eq.(2.24)  with  respect  to  x  to  give 

(2.25)  £x  —  C/xxi 

where  C  is  a  constant  of  the  integration.  Equation  (2.25)  has  one  real  and  two  complex  roots.  Since  we  seek 
only  real  roots  the  complex  roots  are  not  considered.  Taking  into  account  the  boundary  conditions  Eq.(2.2) 
the  above  equation  can  readily  be  integrated,  that  gives 

j(fxx)1/3dx 

(2.26)  e(*)  =  7 - • 

${fxx)lt3dx 

a 

If  a  grid  is  generated  in  accordance  with  the  optimal  mapping  Eq.(2.26)  the  leading  term  of  the  truncation 
error  is  zero  for  all  points  in  [a,  b ]  and  the  global  order  of  accuracy  is  increased  from  2  to  3. 

The  optimal  grid  point  distribution  defined  by  Eq.(2.26)  can  be  applied  if  fxx  is  a  positive  function 
otherwise  the  mapping  becomes  singular  that  leads  to  the  grid  degeneration.  However,  this  problem  can  be 
overcome.  For  that  purpose  we  divide  the  interval  [a,  b]  on  subintervals  where  fxx  is  of  constant  signs.  Let 
x\  <  x  <  X2  be  an  interval  where  the  second  derivative  is  negative,  i.e.  fxx  =  -|/xx|  <  0-  Then,  Eq.(2.26) 
becomes 

f(fxx)1/3dx  -  f  \fxx\1/3dx  f  I fxx\^3dx 

X.  Xi  ®i 

(2.27)  «*)  -  jz - =  — - =  - 

jf  (fxx)1/3dx  -  f  \fBX\V*dx  f  \fxx\1/3dx 
x~  *1  •  X1 

From  Eq.(2.27)  it  follows  that  the  metric  coefficient  fx  given  by  Eq(2.26)  is  strictly  positive  in  the  interval 
where  fxx  is  negative.  Taking  into  account  the  fact  that  the  same  formula  Eq.(2.27)  remains  valid  for  positive 
fxx  the  intervals  of  positive  and  negative  signs  except  for  the  inflection  points  of  the  function  f(x)  can  be 
joined  so  that 

E  f  \fxx\1/3dx 

(2.28)  £(*)  = 

E  /  \fxx\1/3dx 

3  Xj+ 0 

where  Xj  are  the  inflection  points  of  f(x).  To  add  the  inflection  points  fxx(%j )  —  0  to  the  above  integrals 
special  consideration  is  required. 

Let  xq  be  a  point  of  inflection  of  the  function  f(x ),  i.e.  fxx(x 0)  =  0.  Note  in  passing  that  if  we  modify 
the  function  f(x)  by  adding  an  arbitrary  linear  function  the  optimal  grid  Eq.(2.26)  remains  unchanged. 
Furthermore,  if  the  function  f(x)  is  linear  in  the  whole  interval  [a,  b]  then  from  Eq.(2.25)  it  follows  that 
=  0,  Vz  6  [a,  b\.  It  results  in  that  the  grid  step  size  in  the  physical  domain  Ax  =  A£/£x  tends  to  infinity. 
It  can  be  interpreted  as  to  approximate  the  first  derivative  of  the  linear  function  exactly  an  arbitrary  large 
grid  spacing  can  be  used.  Expanding  fxx  in  a  Taylor  series  about  x  =  xQ  in  Eq.(2.26)  and  assuming  that 
fxxxipt 0)  7^  0  yield 

fxx{%)  —  fxxx(% o){x  ~  ^0)  H-  0{(x  —  Xq)  ) 


7 


Substituting  the  above  expression  in  Eq.(2.26)  and  neglecting  both  0((x-x 0)2)  and  higher  order  terms  give 
(2.29)  £,x  =  Cfxxx  (*o)(®  —  xq) 


Letting  x  — >  xo  we  have 


£x(x0)  =  lim  (Cfx 

X-+XQ 


'(x0)(x-x0))1/3  =0 


As  noted  above,  this  kind  of  grid  degeneration  when  the  metric  coefficient  £x  vanishes  does  not  impose  any 
restriction  on  the  grid  step  size  at  the  inflection  point.  Therefore,  in  the  vicinity  of  the  inflection  point  the 
original  second  derivative  fxx  can  be  modified  as 


(2.30) 


f  xx  (x)  — 


\fxx\l 

(/..)2+«2 

2e  ’ 


\fxx\  >  e 
|/xx|  ^  ^ 


where  e  is  a  small  positive  parameter.  From  the  above  consideration  it  follows  that  for  an  arbitrary  /  € 
C2[a,6]  the  optimal  mapping  minimizing  the  leading  truncation  error  term  globally  is 


f(fxx)1/3dx 

(2.31)  fr)  =  T1 - 

J{fxx)1/3dx 


To  estimate  the  asymptotic  truncation  error  of  the  second-order  difference  expression  for  fx  on  the 
optimal  grid  Eq.(2.26)  we  rewrite  Eq.(2.8)  including  the  third-order  terms 


(2.32) 


Lh(fx) 


ft  +  g2A$2/^  +  CsAff^  +  0(A£4) 
*{  +  C2  +  C3A?xf 


Linearizing  Eq.(2.32)  and  collecting  the  terms  of  0( A£2)  and  0(A£3)  the  first  two  leading  terms  in  the 
truncation  error  are 

(2.33)  Ta(aO  =  Ca^f-  [/««*{  -  *«««/«]  +  [/«(4)*«  ~  *«4)/f]  • 

xt  xi 

Since  the  first  term  on  the  right  hand  side  of  Eq.(2.33)  is  vanished  on  the  optimal  grid  defined  by  Eq.(2.26) 
the  asymptotic  truncation  error  becomes 

(2.34)  T2(x)  =  [f^xt  ~  x{4)h]  ■ 

To  determine  the  expression  in  the  square  brackets  we  differentiate  Eq.(2.14)  written  for  p  =  2  with  respect 
to  Thus, 


(2.35) 


^xt  ft££x£t  x\^  ft  ~~ 


Resolving  Eq.(2.14)  with  respect  to  and  substituting  it  in  Eq.(2.35)  give 
(2.36)  j/f  “  x £  ^ ft  =  x£ttx%fx x- 

Using  Eq.(2.36)  the  leading  truncation  error  term  on  the  optimal  grid  Eq.(2.26)  can  be  recast  as 


(2.37) 


T2 (3?)  —  C3 X££tfxx 
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Taking  into  account  the  fact  that  for  the  optimal  grid  Eq.(2.26)  holds,  xss(  can  be  represented  in  terms  of 
the  function  f(x )  and  its  derivatives  as  follows 


(2.38) 


«€«  = 


3£xx  £xxx£x 


^fxxx  3fx  ^  fx 

9C*pxx 


where  C  is  the  integration  constant  in  Eq.(2.25).  Substituting  Eq.(2.38)  into  Eq.(2.37)  the  leading  truncation 
error  term  is  given  by 

/n  exr\\  rp  _  /"*  A  <**3  ^ f XXX  f XX 

(2.39)  T2{x)  =  C3  A£  - tyjzp - ‘ 

This  formula  is  valid  for  all  points  from  the  interval  [a,  b }  except  for  the  inflection  points  of  the  function  f(x). 

Let  us  estimate  the  leading  term  of  the  truncation  error  at  a  point  of  inflection  x0  :  fxx(x 0)  -  0.  Since 
we  have  modified  the  second  derivative  fxx  in  the  neighborhood  of  the  inflection  point  Eq.(2.30)  the  second- 
order  term  in  the  truncation  error  does  not  vanish.  Substituting  Eq.(2.30)  into  Eq.(2.33)  and  neglecting 
higher  order  terms  we  have 


T2(x  0)  =  C2A( 


.2  -/,x/x„(/^)-2/3  +  (/x,)1/3/x 


Letting  x  -4  xo  yields 

(2.40)  '  T2(x 0)  =  C2Ae^^- 

Equation  (2.40)  shows  that,  locally,  near  the  inflection  point  only  the  second  order  of  approximation  can  be 
obtained  on  the  optimal  grid.  Note  that  it  is  not  the  case  if  the  function  f(x)  is  linear  because  then,  any 
second-order  accurate  approximation  of  and  x%  in  Eq.(2.5)  on  an  arbitrary  nonuniform  mesh  gives  us  the 
exact  value  of  fx .  By  virtue  of  the  fact  that  the  number  of  the  inflection  points  is  finite  the  L2  norm  of  the 
second-order  accurate  approximation  of  fx  on  the  optimal  grid  should  provide  super  convergent  results. 

In  regions  where  the  function  f(x)  is  discontinuous  the  above  reasoning  is  not  valid  since  the  first  and 
higher  derivatives  do  not  exist  there.  In  contrast  to  the  inflection  point  in  the  vicinity  of  local  extrema  of 
f(x),  where  fxx  achieves  its  maximum  value,  the  fraction  in  Eq.(2.39)  becomes  very  small  so  that  locally, 
even  a  higher  order  of  accuracy  may  be  obtained. 

Remark  2.1  It  can  readily  be  checked  that  standard  grid  adaptation  criteria  such  as  the  arc  length  of 
the  function  f(x)  and  the  second  derivative  fxx  do  not  globally  minimize  the  leading  term  of  the  truncation 
error.  Actually,  using  the  arc  length  grid  adaptation  criterion  the  following  grid  point  distribution  is  obtained 

f  V1  +  fidx 

(2.41)  m  =  -b 

IV 1  + 


Substituting  Eq.(2.41)  into  Eq.(2.33)  yields 

C  a  ~3fxfxx  +  (1  d~  fx)fxxx 

(2.42)  T2(x)  =  C2A£  - (1  +  f2)2 - * 

The  comparison  of  Eq.(2.42)  with  the  leading  term  of  the  truncation  error  obtained  on  a  uniform  grid,  which 
is 


!?"(»)  =  C2Aefxxx, 
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shows  that  this  grid  point  distribution  may  improve  the  accuracy  locally  near  steep  gradients  of  the  function 
f(x).  At  the  same  time,  in  regions  where  fxx  is  much  grater  than  fx ,  e.g.  near  local  extrema  of  }{x ),  the 
actual  order  of  approximation  may  deteriorate  to  one  or  even  be  worth. 

If  instead  of  the  arc  length  adaptation  criterion  one  redistributes  grid  points  in  accordance  with  the 
second  derivative  fxx  the  leading  term  of  the  truncation  error  is 

(2.43)  Ta(x)  =  -C2A£2^£. 

J  XX 

As  it  follows  from  the  above  formula  in  regions  where  \fxx\  <  y/2  the  local  truncation  error  Eq.(2.43)  is 
always  grater  than  the  asymptotic  truncation  error  on  a  uniform  grid. 

Summarizing  what  has  been  said  above  the  following  conclusions  can  be  drawn.  On  the  one  hand,  the 
standard  grid  adaptation  criteria  do  not  provide  the  superconvergence.  On  the  other  hand,  although,  the 
standard  grid  adaptation  techniques  may  locally  improve  the  accuracy  of  calculation  the  global  truncation 
error  may  become  even  larger  than  that  obtained  on  the  corresponding  uniform  mesh.  Despite  the  fact  that 
the  above  consideration  has  been  performed  for  the  second-order  discretization  the  same  conclusion  can  be 
done  for  higher  order  schemes. 

Remark  2.2  We  shall  now  briefly  describe  an  alternative  way  of  the  solution  of  Eq.(2.21).  Integrating 
Eq.(2.21)  by  parts  and  neglecting  the  0( A£)  term  on  the  right  hand  side  yield 

(2.44)  -  hxti  =C< 

where  C  is  a  constant  of  the  integration.  The  above  equation  is  closed  by  using  the  boundary  conditions 
Eq.(2.2). 

In  order  to  find  the  unknown  constant  C  we  rewrite  Eq.(2.44)  in  the  following  form 


(2.45) 


=  C, 


Taking  into  account  the  fact  that 


fxx  — 


i_d_ 
H  dZ 


Eq.(2.45)  is  reduced  to  Eq.(2.25)  and  the  constant  C  can  easily  be  determined,  that  gives 


(2.46)  C  =  ( J(fxx)l/3dx 

The  boundary  value  problem  Eq.(2.44),(2.46),(2.2)  should  be  solved  numerically.  If  at  some  point  and 
are  equal  to  zero  simultaneously  Eq. (2. 44), (2. 46)  is  degenerated.  The  problem  can  be  overcome  by  modifying 
the  derivatives  and  the  constant  C  as  follows 


f  —  La  = 

*  ~  u  “  (cfxxy/3 


hc  = 


fxxfjx  tjxxfx  _  3{fxx)  fxxxfx 

&  "  C2/3(/as)5/3 
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c=  (j  {fxx)1/3dx 

where  fxx  is  given  by  Eq.(2.30),  fx  and  fxxx  are  calculated  by  differentiating  and  integrating  fxx  with  respect 
to  x ,  accordingly.  Since  the  function  fxx  is  strictly  positive  in  the  entire  computational  domain  the  first 
derivative  /f  is  a  positive  function  as  well.  It  makes  the  modified  equation  fully  consistent  with  Eq.(2.31). 

It  should  be  stressed  that  there  are  several  differential  forms  of  the  optimal  grid  generation  equation. 
For  example,  instead  of  integration  of  Eq.(2.21)  by  parts  we  may  consider  Eq.(2.23)  as  a  differential  equation 
for  the  optimal  grid  point  distribution.  Since  each  of  the  differential  equations  has  its  advantages  and 
disadvantages  at  the  present  time,  it  is  difficult  to  say  which  one  of  them  is  better. 

2.3.  High-Order  Approximations,  p  >  3.  If  /$  and  are  approximated  identically  by  a  third-order 
accurate  formula  the  optimal  grid  generation  equation  written  in  operator  form  in  the  physical  space  is 


(2.47) 


l_cT|4 

Ax  dx~ 


f~fx 


Ax  ®x. 


x  —  o. 


Performing  the  indicated  differentiation  we  have 

(2.48)  fxx  (15&  —  ^(,xf>xxx)  +  (,X  (j~6£xxfxxx  +  £xfX  ^  =  0 

Although,  the  above  equation  is  much  more  complicated  than  the  analogous  one  derived  for  the  second-order 
discretizations  Eq.(2.24)  we  shall  construct  the  solution  of  Eq.(2.48)  in  a  similar  form.  On  the  one  hand, 
a  solution  in  the  form  of  £  ==  g{fx),  where  g  is  an  arbitrary  function  of  fx ,  is  not  appropriate  since  in  this 
case  the  /i4)  term  in  Eq.(2.48)  can  not  be  canceled.  On  the  other  hand,  if  a  solution  depends  on  fxxx  or 
higher  derivatives  of  f(x)  it  gives  rise  to  the  f ^  term  in  the  truncation  error  which  is  not  canceled  as  well. 
Therefore,  we  shall  seek  the  solution  of  Eq.(2.48)  in  a  form  similar  to  Eq.(2.25) 


(2.49) 


e.  =  c(fxxr. 


Substituting  Eq.(2.49)  into  Eq.(2.48)  the  leading  truncation  error  term  can  be  written  as 

(2.50)  T3(x)  =  [a{2  -  11  a){fxxx)2  +  (4a  -  1  )/**/i4)] 

In  contrast  to  the  second-order  discretization,  for  the  third-order  approximation  the  leading  term  of  the 
truncation  error  does  not  vanish  at  any  a  —  const .  Assuming  that  the  parameter  a(x)  is  a  function  which 
weakly  depends  on  x  and  setting  the  leading  truncation  error  term  equal  to  zero  the  following  quadratic 
equation  for  a(x)  is  obtained 

(2.51)  a(x)(2  -  lla(x))(fXxx)2  +  (4a(a;)  -  1  )/x*/|4)  =  0 
The  solution  of  Eq.(2.51)  is 

(2.52) 
with 


ai,2  =  ^j-  (l  +  2r(x )  ±  a/1  -  7 r{x)  +  4 r(x)2)  , 


r(x)  — 


fxxfi 


(4) 


{fxxx)2 
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Without  loss  of  generality  it  is  assumed  that  fxxx  /  0.  If  fxxx  —  0  then  the  solution  of  Eq.(2.51)  is  a  =  1/4. 
Note  that  the  function  a(x )  should  be  positive  in  the  entire  physical  domain  otherwise,  the  mapping  Eq.(2.49) 
with  a  <  0  concentrates  grid  points  where  / (x)  is  linear  and  makes  the  grid  very  coarse  where  the  second 
derivative  fxx  is  large.  Since  the  above  analysis  is  valid  if  the  function  a(r)  slightly  depends  on  x  we  construct 
a  as  follows 


(2.53) 


a(r)  =  / 


yj-  (l  +  2r  +  Vl  -  7r  +  4r2) , 

_  _48_3  ,  18r2  _  JLr  4.  JL. 
343'  ^49'  22'  ~  11’ 

i  (1  +  2r  -  VI  -  7r  +  4r2) , 


r  <  0 
0<r  <  \ 

r>l 


where  the  polynomial  in  Eq.(2.53)  has  been  chosen  so  that  a(r)  is  a  continuously  differentiable  function  of 
r.  A  plot  of  a  versus  r  is  shown  in  Fig.2.1.  As  can  be  seen  in  the  figure,  the  function  a(r)  is  practically 
equal  to  1/4  in  the  whole  range  of  r  except  for  an  interval  -1  <  r  <  3.  Although,  a(r)  is  quite  smooth  the 
function  a(x)  may  be  non-smooth  because  it  depends  on  fxx,  fxxx  and  fx  ^  which  are  calculated  numerically 
and  may  therefore  be  very  oscillatory.  In  numerical  applications  the  function  a(x)  should  be  smoothed  to 
meet  the  requirements  used  for  the  derivation  of  Eq.(2.50). 


Fig.  2.1.  Parameter  a  for  a  third-order  accurate  discretization. 


Such  a  choice  of  a(x)  provides  that  the  leading  truncation  error  term  is  approximately  equal  to  zero 
in  the  entire  physical  domain.  As  it  follows  from  Eq.(2.49),  the  second  derivative  fxx  must  be  a  positive 
function  on  [a,  b\.  Note  that  a  general  property  of  both  Eq.(2.47)  and  Eq.(2.14)  is  that  if  is  a  solution 
of  Eq.(2.47)  then  is  a  solution  of  Eq.(2.47)  as  well.  The  same  is  true  for  the  function  f(x)  and  its 
derivatives ,  i.e.  if  we  substitute  fxx  =  —  fxx  h^0  Eq.(2.47)  we  get  the  same  equation  in  terms  of  fxx.  Hence, 
the  second  derivative  fxx  in  Eq.(2.49)  can  be  replaced  with  Eq.(2.30).  Thus,  if  /f  and  x £  are  evaluated  by 
the  same  third-order  accurate  formula  the  optimal  grid  point  distribution,  which  minimizes  the  leading  term 
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of  the  truncation  error  in  the  entire  computational  domain,  is 


J(fxxja{x)dx 

(2.54)  £(x)  =  -b - , 

f{fxx)a{x)dx 

a 

where  fxx  and  a(x)  are  defined  by  Eq.(2.30)  and  Eq.(2.53),  respectively. 

From  the  above  analysis  one  can  see  that  the  same  strategy  used  for  the  third-order  approximation  can 
be  applied  to  higher  order  discretizations.  Actually,  the  leading  term  of  the  truncation  error  for  an  arbitrary 
pth-order  approximation  of  fx  is 

(2.55)  T,(0  = 

Accounting  for  the  following  relations  between  the  f-  and  ^-derivatives  written  in  operator  form 

d_  _  }_d_ 
d(,  ~  ix  dx 

dn  _ 

d£n~%txdx. 

the  truncation  error  can  be  transformed  into  the  physical  space  as  follows 


Tp(x)  =  Cp  Aetx 


Expanding  the  power  of  the  derivative  operator 


f-fx 


(2  57,  [«r /-[[**]*]- 

-  [**]'”  -«+<»+ 1)  [**]'-  [**]  *£  +  - + [t*]  -  [**]’ « 

it  can  be  seen  that  the  term  with  fx  in  Eq.(2.56)  is  canceled  and  therefore  the  highest  derivatives  of  £(#) 
and  f(x)  in  the  truncation  error  Tp(x)  are  and  /iP+1),  respectively.  Assuming  that  on  the  optimal  grid 
the  leading  term  of  the  truncation  error  is  of  the  order  of  0(A£)  we  shall  seek  £(x)  as  a  function  of  f(x) 
and  its  derivatives.  Comparing  the  highest  derivatives  of  £  and  /  one  can  observe  that  if  £x  —  g{f,fx) 
then  the  term  /iP+1^  in  Eq.(2.57)  is  never  canceled  while  if  is  a  function  of  fxn\  n  >  3  it  introduces  the 
uncancellable  /in+p_1)  term  in  the  truncation  error  Tp(x).  In  a  similar  manner  to  the  first-,  second-,  and 
third-order  approximations  the  optimal  grid  for  the  pth-order  accurate  discretization  is  sought  in  the  form 
of  Eq.(2.49).  Substituting  Eq.(2.49)  into  Eq.(2.57)  the  leading  truncation  error  term  becomes 

(2.58)  Tp(x)  =  ^[1  —  +  1)]  fiP+l ^  +  ocG(a ,  fxx,  fxxx,  *  •  •  > 

In  the  above  formula  it  has  already  been  taken  into  account  that  the  second  term  on  the  right  hand  side  is 
proportional  to  a.  This  is  no  surprise  since  for  a  =  0,  which  corresponds  to  a  uniform  mesh,  the  asymptotic 
truncation  error  Tp(x )  is  reduced  to  CpA^/iP+1)  that  is  why  all  the  terms  in  Eq.(2.58)  except  for  /iP+1) 
have  to  be  proportional  to  ct.  For  example,  for  fourth-  and  fifth-order  discretizations  the  leading  truncation 
error  terms  obtained  on  the  optimal  grid  Eq.(2.49)  are 

(2.59)  T4(x)  =  ((1  -  5 a) /<5>  +  a  j-10a(l  +  5 aO^rJr  +  5(9a  " 
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and 


Tt(x)  =  ((1  -  6a) /i6)  +  a  {(6  +  49a  +  196a2  +  274a3)g^- 

^'U)>  -(7  + 97a  +  421aa)(/f^Vjf-  +  3(27a  -  2)^^  +  2(26a  -  , 


respectively.  As  it  follows  from  Eq.(2.58)  at  any  a  =  const  both  terms  on  the  right  hand  side  do  not  vanish 
simultaneously.  To  minimize  the  leading  term  of  the  truncation  error  the  following  procedure  is  proposed. 
At  each  grid  point  the  parameter  a  is  found  as  the  solution  of  the  nonlinear  equation  T(a)  =  0,  which  is 
solved  by  the  Newton’s  method.  That  choice  of  a  provides  that  the  leading  truncation  error  term  is  vanished 
on  the  optimal  grid.  Since  the  above  consideration  is  valid  only  if  a  slightly  depends  on  x  the  function  a(x) 
has  to  be  smoothed  in  numerical  applications. 

Remark  2.3  If  p  ->-  -boo,  i.e.  the  order  of  approximation  is  infinitely  large  the  leading  term  of  the 
truncation  error  Eq.(2.58)  is  vanished  for  a  -b  0.  In  other  words,  the  higher  is  the  order  of  approximation 
used  to  evaluate  /$  and  the  more  uniform  is  the  grid  which  minimizes  the  leading  truncation  error  term. 
In  the  limit  of  infinitely  high-order  approximations  a  uniform  grid  is  optimal  in  the  sense  of  minimization  of 
the  asymptotic  truncation  error. 


3.  Grid  Adaptation  in  Multiple  Dimensions.  The  two-dimensional  transformation  of  the  first 
derivative  is  given  by 


(3.1) 


fx  = 


Vnh  -  ViU 
J 


where  the  Jacobian  of  the  transformation  is 


J  —  x^y^  x^y^. 

Approximating  the  £  and  7]  derivatives  in  Eq.(3.1)  by  some  pth-  and  gth-order  finite  difference  formulas, 
respectively,  we  get 

_  . ,  ,  (ft,  +  CqA^+1))(ft  +  H))  ~  (j/£  +  CpA^y(f+1)){fv  +  C^Arff^) 

hUx  (*<  +  CpA^x(^+1))(yv  +  C,Aij*»i9+1))  -  (*„  +  C,ArM,+1))(»€  +  CrA?y{*+1)) 

(3.2) 

In  the  above  expression  it  has  already  been  taken  into  account  that  the  metric  coefficients  x^^y^  and  x^^y^ 
are  evaluated  by  the  same  finite  difference  operators  which  are  used  for  calculating  and  /^,  respectively. 

In  view  of  the  fact  that  the  mapping  used  is  nonsingular  J  >  0,  the  denominator  of  Eq.(3.2)  can  be 
linearized  that  yields 

Lh(fx)  =  i  [ft,/{  -  ViU  +  CpA?(ynf(p+1)  -  y{P+1)fv)  +  CtArj“ (fiyi9+1)  -  %49+1))]  x 

(3-3)  jj  _  c^^P+1)  _  yMXti)  _  _  „€a.0ri-i))]  +  0(A?+\  A^1) 

Multiplying  out  the  terms  in  the  square  brackets  and  neglecting  higher  order  terms  the  leading  term  of 
truncation  error  becomes 

rM«,u)  =  $  {cP a?  [v4P+1)  -  y?+1)fv  ~  f*(vAP+i)  ~  + 

(3'4)  +CqArjq  [/^+1)  -  ViA9+1)  -  U(xdq+1)  -  yAv+1))] } 

As  in  the  ID  case,  the  truncation  error  Tp,g  consists  of  two  different  parts,  one  of  which  arises  from  the 
evaluation  of  the  metric  coefficients  x^^y^^xv^yv  and  the  second  one  occurs  due  to  the  approximation  of 
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and  fv.  From  Eq.(3.4)  it  follows  that  if  the  absolute  value  of  the  first  expression  in  the  square  brackets  is  less 
than  0(A£)  and  the  abolute  value  of  the  second  one  is  less  than  O(Arj)  then  the  global  truncation  error  is 
0(A£P+1 ,  Ar)q+1)  rather  than  0(A£P,  At?9).  Thus,  to  increase  the  order  of  the  finite  difference  approximation 
Eq.(3.2)  by  one  globally  grid  points  should  be  redistributed  so  that  the  following  equations  hold 


J  -  a^h)  =  la,  ft  -  Vtf,)  -  *nVtr))  + 

J  -  la, ft  -  vM  (*s»S,+‘>  -  %4'+1))  +  Ol&v) 


Removing  the  parentheses 
(3.6) 


and  rearranging 

y,  |/t’+11 

at  fi'*') 


the  corresponding  terms  Eq.(3.5) 

-  fyV(l+l)  -  fAP+l)]  =  0^J 

-  fyyiq+1)  -  fJv+1)\  =0(AV)J 


can  be  reduced  to 


Note  that  a  reduction  of  the  0( A£)  and  0(Ar/)  terms  in' Eq.  (3.6)  decreases  the  truncation  error  on  the 
optimal  grid. 

The  above  equations  can  be  treated  as  the  optimal  grid  generation  equations  in  the  sense  of  minimization 
of  the  leading  truncation  error  term.  It  should  be  noted,  that  if  2/^  =  0  in  the  entire  computational  domain 
Eq.(3.6)  is  reduced  to  Eq.(2.14).  At  the  same  time,  if  the  y  coordinate  does  not  depend  on  7?,  i.e.  y  =  y(0 
Eq.(3.6)  is  simplified  to 


(3.7) 


xvf(g+1)  -  =  0(Arj)xl, 


that  can  be  treated  as  an  analog  of  Eq.(2.13)  in  the  r)  coordinate. 

Another  very  useful  property  of  the  optimal  mapping  is  that  Eq.(3.6)  are  invariant  with  respect  to  both 
translation  and  stretching  of  the  x,  y  and  £,  rj  coordinates.  Summarizing  the  above  properties  of  Eq.(3.6)  one 
may  conclude  that  the  2D  optimal  grid  generation  equations  are  fully  consistent  with  the  ID  counterpart 
Eq.(2.14). 

The  present  approach  can  directly  be  extended  to  three  dimensions.  Actually,  the  three-dimensional 
transformation  of  the  first  derivative  is 

(3.8)  +  + 

where  the  Jacobian  of  the  mapping  is  given  by 


J  =  xtVrjZt  +  x^zt  +  xcy^zv  -  x^z^  -  x^zt  -  X^ynz^, 


With  pth-,  gth-,  and  rth-order  finite  difference  approximations  for  the  77- ,  and  (-derivatives,  respectively, 
we  have 

/f  x  (^c z^y  ~  ^ySr]z)5^f  -4-  (S^yS^z  -  S^zS^S^f  +  (S^zS^y  -  d^yd^S^f  A£p+1,  A7jq+17  A(r+1), 

h  x  S^xS^yS^z  +  SfjxS^yS^z  d^xS^yS^z  —  S^xS^yS^z  —  S^xS^yS^z  —  S^xS^yS^z 

(3.9) 

where  the  differential  operators  ,  SV)  and  6$  are  defined  by 

h  =  -§i  +  cv^pw±i 

(3.10)  Sri  =  Tftj  +  CqArjq  0)?,+  | 

6(  =  fi+CrAC-§£r. 


Here,  Cv\  Cq ,  and  Cr  are  constants  dependent  on  particular  pth-,  gth-,  and  rth-order  finite  difference 
approximations  which  are  applied  to  discretize  the  77-,  and  (-derivatives,  accordingly.  In  Eq. (3.9), (3. 10) 
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it  has  already  been  accounted  for  that  the  metric  coefficients  are  approximated  by  the  same  finite  difference 
expressions  which  are  used  for  evaluating  /,,,  and  fcr 

Having  linearized  the  fraction  in  Eq.(3.9)  the  leading  truncation  error  term  is  written  as 

aWK.u, 0=7  {cPHv  [^(p+1)  -  %p+1)k]  +  cq^nq  [4?+1)  -  49+1)/*]  + 

(3'n)  +Cr  AO  -  4r+1)/,]  }  , 

where 

Fe(p+1)  =  f^+1)(z(yv  -  ycz„)  +  yf+1\znh  ~  ZC  k)  +  ^+l){ycfr,  ~  Vvk) 

49+1)  =  4,+1W  -  Hvd  +  yrl\nh  -  nk)  +  4,+1W<  -  y<k) 

Ft Jr+1)  =  f\r+l){zr,Vi  -  VvH)  +  4r+1  Wo  -  znk)  +  4r+1W*  “ 

(3-10  jW>  =  ^  +  yf+1\zvX(  -  zcxv)  +  zf+1\ycxn  -  yvx c) 

45+1)  =  4?+ l)(zt2/C  -  +  VvHz Cxi  ~  zHxc)  +  zn+1Hyt;x C  ~  V(xi) 

4r+1)  =  x£+1)(zr,ys  -  ynz<:)  +  y^+l) {z(.xn  -  z^)  +  z^Hy^  -  y(xv) 

Similarly  to  the  ID  and  2D  cases  described  above  the  leading  term  of  the  truncation  error  Eq.(3.11)  can  be 
divided  into  two  parts.  The  first  part,  which  also  exists  on  a  uniform  mesh,  is  due  to  the  approximation  of 
A 5  fr\ ?  and  /^.  The  second  part,  which  is  vanished  on  a  uniform  Cartesian  mesh  is  caused  by  the  evaluation 
of  the  metric  coefficients.  Prom  Eq.(3.11)  it  is  apparent  that  if  a  grid  is  constructed  so  that  the  first  term 
in  the  square  brackets  is  of  the  order  of  0(A£),  the  second  one  is  of  the  order  of  0(A?7),  and  the  third  one 
is  of  the  order  of  0(AQ  for  all  £  €  [0, 1],  r\  G  [0, 1],  and  C  e  [0, 1]  then  the  global  order  of  approximation 
of  the  difference  operator  Eq.(3.9)  in  £,  77,  and  (  on  the  optimal  grid  is  increased  from  p,  q,  and  r  to  p  +  1, 
q  4.  l,  and  r  +  1,  respectively.  Hence,  in  the  sense  of  minimization  of  the  leading  truncation  error  term  the 
grid  adaptation  criteria  are 


(3.13) 

ff+1)  -  SJ\P+l)  =  O(A0  J 

(3.14) 

F^+1)  -  =  0(At])J 

(3.15) 

Fc(r+1)  -  fj£+1)  =  0(A<)J. 

Note  that  the  above  equations  are  not  a  system  of  equations  and  can  be  considered  separately.  If  it  is 
necessary  to  improve  the  accuracy  with  respect  to  the  £  coordinate  alone  a  grid  should  be  generated  so  that 
only  Eq.(3.13)  holds.  However,  if  it  is  desirable  to  increase  the  order  of  approximation  of  fx  by  one  in  the 
£,  77,  and  £  coordinates  simultaneously  then  the  grid  has  to  obey  the  system  of  equations  Eq.(3.13)-(3.15). 

As  in  the  case  of  two  dimensions  the  3D  grid  adaptation  criteria  Eq.(3.13)-(3.15)  can  be  simplified.  After 
the  substitution  of  Eq.(3.12)  in  Eq.(3.13)-(3.15)  and  considerable  algebraic  manipulation  the  grid  adaptation 
equations  can  be  rewritten  in  a  very  compact  form 

(zcv,  -  y(zv)  [/<p+1)  -  Ux\v+l)  -  kVi+1)  -  kzt+1)]  =  °^J 
(3.16)  (ycz€  -  zm)  [fvq+1)  ~  h  49+1)  -  fyVv+1)  ~  fA9+1)\  =  0(Ar,)J 

(zvVi  -  ynz0  [/((r+1)  -  k4+1)  -  fyy[r+l)  -  fAr+1)\  =  °^0J, 

where  /*,  fy,  and  fz  are  the  first  derivatives  with  respect  to  the  x ,  y,  and  0  coordinates,  respectively.  One  of 
the  characteristic  features  of  the  above  equations  is  that  they  do  not  depend  on  the  coefficients  Cp,  Cq ,  and 
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Cr.  Consequently,  if  in  each  spatial  direction  the  metric  coefficients  and  the  first  derivatives  of  /(£,??}£)  are 
evaluated  consistently  by  some  hybrid  finite  difference  operators  then  the  grid  adaptation  criteria  Eq.(3.16) 
can  be  applied  in  the  whole  computational  domain  regardless  of  points  where  the  hybrid  scheme  switches 
from  one  approximation  to  another.  A  comparison  of  Eq.(3.16),  Eq.(3.6),  and  Eq.(2.13)  shows  that  the  3D 
grid  adaptation  criteria  Eq.(3.16)  are  reduced  to  Eq.(3.6)  if  =  zv  =  0,  z^  ^  0,  while  if  in  addition  to  these 
conditions  we  require  that  =  y<;  —  Q,  Vr)  ^  0  Eq.(3.16)  are  reduced  to  the  ID  optimal  grid  generation 
equation  Eq.(2.13).  In  a  similar  manner  as  Eq.(2.13)  and  Eq.(3.6),  it  is  easy  to  prove  that  Eq.(3.16)  are 
invariant  with  respect  to  stretching  and  translation  of  both  the  physical  and  computational  coordinates. 

As  it  follows  from  the  analysis  presented  in  the  foregoing  section  the  grid  adaptation  equation  does  not 
assure  that  the  coordinate  mapping  obtained  as  the  solution  of  Eq.(2.14)  is  not  singular.  Since  Eq.(3.16) 
is  converted  to  Eq.(3.6)  and  in  its  turn  Eq.(3.6)  is  reduced  to  Eq.(2.14)  if  the  dimension  of  the  space  is 
decreased  by  one,  the  same  singularity  may  occur  in  two  and  three  dimensions  as  well. 

Equations  (3.6)  and  (3.16)  have  to  be  closed  by  corresponding  boundary  conditions.  Since  these  equations 
are  (p+  l)th-order  partial  differential  equations  p+ 1  boundary  conditions  should  be  imposed  at  each  couple 
of  the  opposite  boundaries  (i.e.  £  —  0  and  £  =  1;  rj  =  0  and  77  =  1;  £  =  0  and  £  =  1)  to  find  the  unique 
solution.  However,  at  each  boundary  we  have  only  one  boundary  condition.  For  example,  in  the  3D  case  in 
the  £  coordinate  we  have 

(3.17)  .  £(x,y,s)  =  0,  £(x,y,2f)  =  1. 

In  other  words  Eq.(3.6)  and  Eq.(3.16)  are  not  closed.  The  situation  becomes  even  more  uncertain  when 
only  one  of  the  grid  adaptation  criteria  is  used.  However,  this  uncertainty  gives  us  additional  degrees  of 
freedom  and  at  the  same  time,  it  is  conceivable  that  there  exists  more  than  one  optimal  grid  satisfying  the 
criteria  Eq.(3.6)  or  Eq.(3.16).  From  this  standpoint  both  Eq.(3.6)  and  Eq.(3.16)  should  be  treated  as  the 
grid  adaptation  criteria  rather  than  the  optimal  grid  generation  equations. 

One  of  the  most  general  structured  grid  generation  strategies  is  based  on  the  variational  approach 
proposed  by  Brackbill  and  Saltzmann  in  [17].  In  this  method  a  grid  is  generated  as  the  solution  of  the 
minimization  problem.  By  forming  the  variational  principle  using  a  linear  combination  of  the  integral 
measures  of  smoothness,  orthogonality,  and  adaptation,  a  system  of  elliptic  equations  is  derived.  The  new 
grid  adaptation  criteria  can  be  incorporated  into  this  approach  by  constructing  an  integral  measure  of 
adaptation  so  that  the  Euler-Lagrange  equations  associated  with  the  minimization  of  this  integral  alone  give 
us  Eq.(3.16).  On  the  one  hand,  the  minimax  principle  guarantees  that  the  coordinate  mapping  obtained  as 
the  solution  of  this  minimization  problem  is  not  singular.  On  the  other  hand  ,  the  new  grid  adaptation  criteria 
provide  that  the  leading  term  of  the  truncation  error  is  minimized  so  that  the  finite  difference  approximation 
Eq.(3.9)  calculated  on  the  optimal  grid  exhibits  superconvergence  properties. 

Remark  3.1  In  spite  on  the  fact  that  the  present  analysis  has  been  performed  for  the  first  derivative 
fx  it  can  be  directly  extended  to  an  equation  or  a  system  of  equations,  which  can  be  represented  as 

(3.18)  fx(x)=d(x). 

For  example,  for  the  steady  state  ID  Burgers  equation  written  in  conservation  law  form  we  have 


where  \x  is  a  positive  constant.  A  comparison  of  Eq.(3.19)  and  fx  shows  that  for  the  Burgers  equation  the 
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optimal  grid  can  be  constructed  using  Eq.(2.54)  with 


(3.20) 


It  should  be  pointed  out  that  the  above  conclusion  is  valid  if  the  second  derivative  uxx  =  (ux)x  and  the 
convective  term  (u2  / 2) x  are  approximated  consistently. 

The  same  approach  can  be  applied  to  the  Euler  and  Navier- Stokes  equations.  The  ID  Euler  and  Navier- 
Stokes  equations  can  be  written  in  conservation  law  form  as 


(321)  &  =  °' 

where  F  is  the  inviscid  flux  F*  for  the  Euler  equations  and  F*  —  Fv,  where  Fv  is  the  viscous  flux,  for  the 
Navier-Stokes  equations.  As  it  follows  from  Eq.(3.16),  any  component  of  the  vector  F  can  be  chosen  as  a 
function  with  respect  of  which  a  grid  is  adapted.  Although,  that  choice  provides  increase  in  accuracy  for 
this  particular  vector  component  but  it  may  not  result  in  decrease  in  the  truncation  error  for  the  remaining 
vector  components.  In  fact,  as  there  are  components  of  the  vector  F  as  many  the  optimal  grids  can  be 
generated.  Since  the  different  vector  components  may  have  strong  gradients  and  local  extrema  in  different 
regions  of  the  physical  domain  this  kind  of  grid  adaptation  is  not  effective.  If  this  is  the  case  the  function 
f(x)  can  be  obtained  by  using  the  method  of  least  squares.  Because  of  the  optimal  grid  generation  equations 
are  invariant  with  respect  to  stretching  of  the  function  f(x)  the  vector  components  Fn ,  xi  =  1  ,N  can  be 
normalized  as 


(3.22) 


Fn(aj)  — 


\Fn(x)\ 
max  \Fn(x)\ 

X 


It  results  in  that  all  of  the  vector  components  are  of  the  same  order  and,  consequently,  make  proportional 
contributions  to  the  function  f(x ).  The  resulting  function  f(x )  is  obtained  as  the  solution  of  the  following 
minimization  problem 


(3.23) 


EE  (*•  {Xi)  -  /(Xi)'j 


-4  min 


i= 0  n— 1 

in  the  least  square  sense.  The  function  /  constructed  in  this  fashion  allows  one  to  generate  a  grid  which  is 
optimal  for  the  whole  vector  F  rather  than  for  its  particular  component.  Note  that  the  power  in  Eq.(3.23) 
should  be  chosen  in  accordance  with  the  power  of  the  Lk  norm  in  which  the  solution  of  the  Euler  or  Navier- 
Stokes  equations  is  sought. 


4.  Results  and  Discussion.  To  validate  the  applicability  and  efficiency  of  the  new  method  several 
ID  and  one  2D  test  examples  are  considered.  For  each  ID  test  function  five  series  of  calculation  on  different 
grids  with  the  same  number  of  grid  points  have  been  executed.  The  first  one  is  done  on  a  uniform  grid. 
The  second  one  uses  the  standard  grid  adaptation  criterion  based  on  the  arc  length  or  the  second  derivative 
of  the  test  function.  The  third  one  is  performed  on  the  optimal  grid  obtained  as  the  analytical  solution  of 
Eq.(2.14).  The  fourth  one  employs  the  optimal  grid  Eq.(2.54)  generated  numerically  by  using  the  following 
approximation  for  the  second  derivative 


(4.1) 


( fxx)i  — 


hjfj+i  —  (hi  +  frj+i)/i  +  fj—i 

1)/2 


X{  Xi— i, 


which  is  reduced  to  the  second-order  three-point  central  approximation  of  fxx  if  an  equispaced  grid  in  the 
physical  domain  is  used.  The  integrals  in  Eq.(2.54)  is  computed  using  the  trapezoidal  rule  integration.  As 
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a  result  of  this  integration  the  strictly  increasing  function  £(x)  is  obtained  which  is  then  reversed  by  using 
a  third-order  accurate  piecewise  spline  interpolation.  The  fifth  calculation  is  also  executed  on  the  uniform 
grid,  however,  instead  of  a  pth-order  approximation  a  (p  +  l)th-order  accurate  dicretization  is  applied  to 
calculate  both  fa  and  x4.  At  each  boundary  one-sided  pth-order  differences  are  used  for  /?  and  x$. 

In  order  to  estimate  the  accuracy  of  the  method  the  pth-order  finite  difference  approximation  of  fx  is 
compared  with  the  exact  value  of  the  first  derivative  calculated  at  the  same  grid  node  in  the  L2  norm.  The 
order  of  approximation  is  estimated  on  successively  refined  grids  the  coarsest  one  of  which  contains  20  cells 
and  the  finest  one  has  2560  cells. 

4.1.  ID  Test  Examples.  Second-order  approximation ,  p  =  2 

The  first  test  example  is  evaluation  of  the  first  derivative  of  f(x)  =  xm,  0  <  x  <  1  by  using  a  second- 
order  central  differences  for  /f  and  x £.  When  m  is  sufficiently  large  this  function  has  a  boundary  layer  of 
width  0(l/m)  near  x  =  1.  For  this  test  case  the  exact  optimal  grid  point  distribution  defined  by  Eq.(2.21) 
can  be  found  analytically,  which  is 

(4.2)  ^(0  =  ^- 

In  contrast  to  [9]  the  new  grid  adaptation  criterion  provides  the  concentration  of  grid  nodes  near  the  boundary 
layer  of  the  function  f(x ). 


-0 -  Uniform,  2nd  order 

-m -  Numerical  optimal 


Fig.  4.1.  Error  convergence  for  a  second-order  approximation  of  fx,  f(x)  =  x7n  calculated  on:  1)  uniform  grid,  2)  optimal 
grid  generated  numerically ,  3)  analytical  optimal  grid}  4)  grid  adapted  in  accordance  with  the  arc  length  criterion,  5)  uniform 
grid  using  third- order  accurate  discretization. 

An  error  convergence  plot  for  this  test  function  is  presented  in  Fig. 4.1.  As  one  might  expect,  the  L<i 
norm  of  the  truncation  error  calculated  on  a  uniform  grid  exhibits  the  0(A£^)  convergence  rate  which  is 
consistent  with  the  second  order  of  accuracy  of  the  central  differences.  However,  the  same  second-order 
approximation  of  fx  on  the  optimal  grid  Eq.(4.2)  exhibits  the  convergence  rate  which  is  even  higher  than 
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0(A£3).  Although,  the  accuracy  of  fx  obtained  on  the  adaptive  grid  Eq.(2.26)  with  fxx  evaluated  by  Eq.(4.1) 
is  slightly  less  compared  to  the  optimal  grid  Eq.(4.2)  results  the  order  of  approximation  is  about  3.5.  To  show 
the  superiority  of  the  present  method  over  the  standard  grid  adaptation  criterion  Eq.(2.41)  the  truncation 
error  calculated  on  grids  adapted  in  accordance  with  the  arc  length  of  f(x)  is  also  shown  in  Fig.4.1.  In 
spite  of  the  fact  that  the  standard  grid  adaptation  technique  slightly  improves  the  accuracy  of  calculation  in 
comparison  with  the  equispaced  grid  point  distribution  the  convergence  rate  is  less  than  0(Af2).  We  want 
to  emphasize  that  the  new  grid  adaptation  criterion  Eq.(2.26)  provides  not  only  superconvergent  results,  but 
on  the  finest  mesh  it  reduces  the  error  by  6  orders  of  magnitude  compared  to  the  uniform  grid  results. 

An  advantage  of  the  consistent  grid  adaptation  Eq.(2.14),  which  is  based  on  the  fact  that  the  truncation 
errors  due  to  the  approximation  of  and  cancel  each  other,  becomes  obvious  when  the  optimal  grid  results 
are  compared  with  those  obtained  by  using  a  third-order  accurate  approximation  on  a  uniform  grid.  Figure 
4.1  shows  that  both  the  second-order  approximation  on  the  optimal  grid  and  the  third-order  discretization 
on  the  uniform  grid  with  the  same  number  of  grid  points  reveal  the  0(A£3)  convergence  rate.  However,  the 
optimal  grid  results  are  about  103  times  more  accurate. 

It  should  be  noted  that  the  optimal  grid  Eq.(4.2)  is  essentially  non-smooth  and  does  not  meet  the 
standard  criterion  of  smoothness,  which  is  \x^/x^\  <  0(1)  [18].  Furthermore,  the  optimal  mapping  Eq.(4.2) 
is  singular  at  the  point  £  =  0  where  x^  —>  oo.  In  spite  on  this  fact,  the  above  comparisons  corroborate  the 
theoretical  analysis  and  demonstrate  the  advantage  of  the  new  grid  adaptation  criterion  over  the  standard 
approaches. 


■© -  Uniform,  hybrid 


FlG.  4.2.  Error  convergence  of  a.  second-order  hybrid  approximation  calculated  with  the  consistent  and  inconsistent  dis¬ 
cretizations  of  the  metric  coefficient  on  the  optimal  and  uniform  grids . 

Another  very  useful  characteristic  feature  of  the  new  method  is  its  generality  in  the  sense  that  if  a 
single  second-order  hybrid  discretization  is  used  for  both  and  x$  the  same  optimal  mapping  Eq.(4.2) 
minimizes  the  leading  truncation  error  term.  To  demonstrate  this  property  the  error  convergence  of  the 
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hybrid  approximation  obtained  on  the  uniform  and  optimal  grids  with  the  same  number  of  grid  points  are 
depicted  in  Fig.4.2.  The  hybrid  difference  operator  is  constructed  as  follows 


u  fW)  =  I  ~  \ even 

l^/j  \  +  ^/i+l  —  /*+2)>  *  °dd 

The  identical  approximation  is  employed  for  the  metric  coefficient  x$ .  A  comparison  shows  that  the  global 
order  of  the  consistent  approximation  of  and  x £  is  increased  by  one  on  the  same  optimal  grid  Eq.(4.2) 
used  for  the  non-hybrid  approximation.  As  has  been  shown  in  Section  2,  the  approximation  of  the  metric 
coefficient  and  the  first  derivative  /$  should  be  the  same  otherwise  the  optimal  mapping  defined  by  Eq.(2.26) 
does  not  minimize  the  leading  truncation  error  term.  To  show  that  the  discretiztion  of  the  metric  coefficient 
plays  a  crucial  role  in  reduction  of  the  truncation  error  we  approximate  x%  by  a  two-point  central  difference 
expression  in  the  whole  computational  domain  and  use  the  same  hybrid  scheme  Eq.(4.3)  for  /^.  An  error 
convergence  plot  for  this  inconsistent  approximation,  which  is  also  depicted  in  Fig.4.2,  shows  that  if  the 
metric  coefficient  are  evaluated  in  a  different  way  than  the  order  of  approximation  on  the  optimal  grid 
deteriorates  to  2  as  well  as  the  truncation  error  increases  by  a  factor  of  103  compared  to  the  consistent 
discretization  results. 

The  second  test  function  considered  is 


(4.4) 


/(*)  = 


1 

(em  -  l)x  +  1 5 


0  <x  <  1. 


-0 -  Uniform,  2nd  order 

-■ -  Numerical  optimal 


Fig.  4.3.  Error  convergence  for  a  second-order  approximation  of  fx,  f(x)  —  l/((e7n  —  l)x  + 1)  calculated  on:  1)  uniform 
grid,  2)  optimal  grid  generated  numerically,  3)  analytical  optimal  grid,  4)  grid  adapted  in  accordance  with  the  arc  length 
criterion,  5)  uniform  grid  using  third-order  accurate  discretization,  6)  numerical  optimal  grid  generated  iteratively. 

In  the  present  test  example  the  parameter  m  was  chosen  to  be  5.  This  function  has  a  boundary  layer 
of  width  0(m/(em  -  1))  at  x  =  0.  For  this  function  the  optimal  grid  generation  equation  Eq.(2.14),  which 
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depends  on  the  order  of  approximation  rather  than  on  a  particular  type  of  discretization,  can  be  solved 
analytically,  that  gives 

-  l 

(4.5)  %opt{£)  ~  em  _  i  ‘ 

It  should  be  emphasized  that  Eq.(2.26)  yields  the  same  optimal  mapping  as  Eq.(4.5).  The  optimal  grid 
Eq.(4.5)  is  the  well-known  exponential  coordinate  transformation,  which  is  widely  used  in  the  literature 
[1],  [18]  for  solving  boundary  layer  problems.  However,  the  mapping  Eq.(4.5)  is  optimal  only  for  a  special 
class  of  functions  such  as  Eq.(4.4)  and  not  optimal  for  other  functions.  Similarly  to  Fig. 4.1  and  4.2,  error 
convergence  plots  for  the  symmetric  second-order  and  hybrid  discretizations  Eq.(4.3)  are  depicted  in  Fig.4.3 
and  4.4,  respectively.  It  is  apparent  in  these  figures  that  the  error  obtained  on  the  optimal  grid  reveals 
the  convergence  rate  of  0(A£3,5)  that  is  even  higher  than  it  follows  from  the  theoretical  analysis.  The 
optimal  grid  point  distribution  constructed  by  the  numerical  integration  of  Eq. (2.26)  reduces  the  truncation 
error  by  about  four  orders  of  magnitude  compared  to  the  uniform  grid  results,  but  it  does  not  provide 
the  same  accuracy  as  the  optimal  grid  Eq.(4.5).  The  accuracy  can  be  improved  if  the  following  iterative 
procedure  is  applied.  Since  the  fxx  approximation  Eq.(4.1)  depends  on  the  grid  spacing  in  the  physical 
domain,  the  second  derivative  can  be  updated  when  the  new  grid  point  distribution  is  found.  For  this  test 
problem  about  15-20  iterations  were  needed  to  reach  the  convergence.  No  attempt  was  made  to  optimize  the 
iteration  process.  Referring  to  Fig.4.3  one  can  see  that  this  procedure  considerably  increases  the  accuracy 
and  provides  practically  the  same  convergence  rate  as  for  the  analytical  optimal  grid  Eq.(4.5). 


-© -  Uniform,  hybrid 


Fig.  4.4.  Error  convergence  of  a  second-order  hybrid  approximation  calculated  with  the  consistent  and  inconsistent  dis¬ 
cretizations  of  the  metric  coefficient  on  the  optimal  and  uniform  grids . 

The  importance  of  the  metric  coefficient  evaluation  is  illustrated  in  Fig. 4. 4.  Analogously  to  the  foregoing 
test  case,  the  inconsistent  discretization  of  and  leads  to  decrease  in  both  the  order  and  accuracy  of 
the  approximation.  When  the  metric  coefficient  and  the  first  derivative  are  evaluated  by  using  the  same 
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hybrid  operator  Eq.(4.3)  the  convergence  rate  obtained  on  the  optimal  grid  Eq.(4.5)  becomes  0( A£3). 

Prom  the  present  theoretical  analysis  it  follows  that  the  new  grid  adaptation  strategy  may  be  quite 
sensitive  to  the  inflection  points  of  the  function  f(x).  In  order  to  verify  this  conclusion  the  following  function 

(4.6)  f(x)  =  [sin(3m.3:)  -  27  sin(ma;)] ,  0  <  x  <  tt, 

which  has  m  inflection  points  has  been  chosen  as  a  test  function.  Despite  the  presence  of  the  inflection  points 
where  fxx  =  0  it  is  possible  to  construct  the  optimal  mapping  analytically  without  using  Eq.(2.30).  It  can 
be  done  if  the  optimal  grid  Eq(2.26)  is  generated  in  each  interval  of  constant  signs  of  fxx  separately.  Thus, 
we  have 

(4.7)  £opt(£)  =  ~0'  —  1)  +  ~  arccos  [2j  —  2m£  -*  1] ,  j 

In  numerical  calculations  the  parameter  m  was  taken  to  be  5.  The  above  optimal  coordinate  transformation 
obeys  Eq.(2.26)  in  the  entire  physical  domain  except  for  the  inflection  points. 


-0 -  Uniform,  2nd  order 

-■ -  Numerical  optimal 


Fig.  4.5.  Error  convergence  for  a  second-order  approximation  of  fx ,  calculated  on:  1)  uniform  grid,  2)  optimal  grid 
generated  numerically,  3)  analytical  optimal  grid,  4)  grid  adapted  in  accordance  with  the  arc  length  criterion,  5)  uniform  grid 
using  third- order  accurate  discretization. 

Figures  4.5  and  4.6  are  analogous  to  Fig.4.1  and  4.2,  accordingly.  As  one  can  see  in  Fig.4.5  the  presence 
of  the  inflection  points  results  in  that  the  convergence  rate  is  0( A£2-5)  that  is  lower  than  it  is  predicted 
from  the  theoretical  analysis.  Nevertheless,  the  optimal  grid  adaptation  reduces  the  truncation  error  by  a 
factor  of  20  compared  to  the  uniform  grid  results.  One  of  the  reasons  of  such  a  behavior  is  the  fact  that 
high-order  derivatives  of  the  function  f{x)  Eq.(4.6)  are  well  bounded  that  makes  the  approximation  of  fx 
on  the  uniform  grid  sufficiently  accurate.  The  use  of  the  standard  grid  adaptation  criterion  based  on  fxx 
Eq.(2.30)  leads  to  deterioration  of  the  convergence  rate  to  0(A£15)  and  at  the  same  time,  the  L2  norm  of 
the  truncation  error  is  about  50  times  less  accurate  than  the  uniform  grid  results.  Figure  4.6  shows  that  the 
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inconsistent  approximation  of  /f  and  x £  increases  the  truncation  error  by  5  orders  of  magnitude  compared 
to  the  consistent  approximation  results  calculated  on  the  optimal  grid. 


-0 -  Uniform,  hybrid 


Fig.  4.6.  Ettot  convergence  of  a  second-order  hybrid  approximation  calculated  with  the  consistent  and  inconsistent  dis¬ 
cretizations  of  the  metric  coefficient  on  the  optimal  and  uniform  grids. 


To  gain  greater  insight  into  where  the  maximum  error  occurs  pointwise  error  distributions  obtained 
on  both  the  uniform  and  optimal  grids  are  shown  in  Fig.4.7.  As  expected,  the  truncation  error  calculated 
on  the  optimal  grid  achieves  its  maximum  values  at  the  inflection  points,  while  the  error  on  the  uniform 
grid  occurs  at  points  where  the  third  derivative  \fXxx\  is  large.  In  contrast  to  the  uniform  grid,  the  most 
accurate  approximation  of  the  first  derivative  fx  on  the  optimal  grid  is  near  the  local  extrema  of  f(x).  For 
demonstrating  the  gain  in  accuracy  in  the  vicinity  of  the  inflection  points  due  to  the  use  of  Eq.(2.30)  instead 
of  fxx  a  pointwise  error  plot  obtained  in  this  case  is  also  presented  in  Fig.4.7.  It  is  significant  that  the  error 
distribution  obtained  on  the  optimal  grid  is  essentially  nonuniform  that  gives  an  indication  of  the  difference 
between  the  present  and  equidistribution  grid  adaptation  criteria. 

From  the  practical  point  of  view  it  is  very  important  to  improve  the  accuracy  of  calculation  if  the  function 
f(x)  is  discontinuous.  In  spite  of  the  fact  that  the  present  analysis  is  not  valid  at  discontinuities  of  f(x)  it 
can  be  used  if  the  discontinuous  function  is  approximated  by  some  smooth  one.  In  this  test  example  the 
following  smooth  function 


(4.8) 


/(*)  = 


■  t17  +  ±  55<">‘  +  15(“gl  +  -  arctan(ac),  -1<«<1 

7 r 


157r(l  +  (ex)2)4 


is  considered  as  a  fitting  of  a  step  function.  In  this  calculation  the  parameter  e  was  taken  to  be  103  that 
results  in  that  the  function  Eq.(4.8)  has  a  pronounced  interior  layer  of  width  0(l/e)  at  x  =  0.  This  function 
has  been  chosen  so  that  the  optimal  grid  point  distribution  Eq.(2.26)  can  be  integrated  analytically.  As  in 
the  foregoing  example,  the  singularity  in  the  optimal  mapping  Eq.(2.26)  due  to  the  inflection  point  at  x  =  0 
can  be  overcome  by  generating  the  optimal  grid  in  the  —  0.5  <  x  <  0  and  0  <  x  <  0.5  intervals  separately, 
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Uniform 


Fig.  4.7.  Pointwise  error  distribution  for  a  second-order  approximation  calculated  on  the  analytical  optimal ,  numerical 
optimal ,  and  corresponding  uniform  grids. 


that  gives 


(4.9) 


xopt(£)  — 


-yf&fv  o<£<0.5 

vtSS T>  0.5<^<l. 


The  error  convergence  of  the  symmetric  second-order  discretization  of  fx  evaluated  on  the  optimal  grid 
Eq.(4.9)  is  compared  with  results  obtained  by  second-  and  third-order  approximations  on  a  uniform  grid  as 
well  as  with  the  truncation  error  calculated  on  grids  generated  by  using  the  standard  Eq.(2.41)  and  new 
Eq. (2. 26) ,(2.30)  grid  adaptation  criteria  in  Fig. 4. 8.  Because  of  the  internal  layer  thickness  is  comparable 
with  the  finest  grid  spacing  none  of  the  uniform  grids  considered  can  provide  second-order  results.  For 
the  analytical  optimal  grid  the  convergence  rate  is  of  the  order  of  0(A£2’5).  Although,  it  is  less  than  the 
theoretical  limit  the  truncation  error  on  the  finest  mesh  (2560  cells)  has  been  reduced  by  more  than  5  orders 
of  magnitude  compared  to  the  uniform  grid  results.  Since  the  standard  grid  adaptation  criterion  Eq.(2.41), 
which  is  widely  used  to  improve  the  resolution  near  steep  gradients  of  the  solution,  does  not  provide  the 
cancellation  of  the  leading  truncation  error  term  these  results  are  about  2  orders  of  magnitude  less  accurate 
than  those  obtained  on  the  optimal  grid  Eq.(2.26), (4.1), (2.30)  as  is  evident  in  Fig.4.8. 

A  comparison  of  the  hybrid  approximation  Eq.(4.3)  on  different  grids  and  using  different  approximations 
for  the  metric  coefficient  x £  is  presented  in  Fig. 4. 9.  If  /f  and  x%  are  evaluated  identically  the  same  optimal 
grid  Eq.(4.9)  provides  superconvergent  results,  while  if  these  approximations  are  different  the  convergence 


rate  is  even  less  than  0(A£2). 


High- order  approximations ,  p  >  3 

For  a  third-order  discretization  the  optimal  grid  generation  equation  Eq.(2.51)  can  not  be  solved  ana¬ 
lytically,  however,  the  solution  can  be  found  in  the  approximate  form  of  Eq.(2.53),(2.54).  The  same  function 
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-Q -  Uniform,  2nd  order 


-■ -  Numerical  optimal 


Fig.  4.8.  Error  convergence  for  a  second-order  approximation  of  fx,  calculated  on:  1)  uniform  grid ,  2)  optimal  grid 
generated  numerically ,  3)  analytical  optimal  grid,  4)  grid  adapted  in  accordance  with  the  arc  length  criterion,  5)  uniform  grid 
using  third-order  accurate  discretization . 


Eq.(4.4),  which  has  been  used  in  the  second  example  is  taken  as  a  test  function.  The  first  derivative  and 
the  metric  coefficient  are  evaluated  by  a  third-order  accurate  formula  as 

(4.10)  (gz)i  =  +  %+ 1  ”  9i+ 2) , 

where  g(£)  is  either  /(£)  or  x(f). 

Figure  4.10  shows  error  convergence  plots  obtained  on  the  optimal  Eq.(2.53),(2.54)  and  uniform  grids 
with  the  same  number  of  grid  cells.  Although,  for  the  mapping  Eq. (2. 53), (2. 54)  the  leading  term  of  the 
truncation  error  is  approximately  equal  to  zero  the  error  convergence  rate  obtained  on  the  optimal  grid  is 
about  0(A£3-8)  that  corroborates  the  theoretical  results.  Note  that  the  same  iterative  technique  used  earlier 
for  the  second-order  approximations  can  be  applied  in  the  present  case  as  well.  However,  due  to  the  fact  that 
the  optimal  coordinate  transformation  Eq.(2.53),(2.54)  is  the  approximate  solution  of  Eq.(2.51)  the  iterations 
do  not  practically  improve  the  accuracy  of  calculation  and  therefore,  these  results  are  not  presented  here. 

The  truncation  error  can  be  reduced  if  the  optimal  grid  generation  equation  Eq.(2.48)  is  solved  numer¬ 
ically.  To  avoid  the  solution  of  the  third-order  differential  equation  a  new  dependent  variable  u(x)  —  £x  is 
introduced.  Then  Eq.(2.48),  which  is  a  second-order  differential  equation  in  terms  of  u(x),  is  integrated  nu¬ 
merically  on  a  uniform  grid  constructed  in  the  physical  domain.  To  close  Eq.(2.48)  the  metric  coefficient  is 
taken  to  be  proportional  to  ( fxx)1/A  at  the  boundaries.  The  metric  coefficient  found  this  way  is  integrated 
and  the  optimal  grid  point  distribution  is  obtained  by  a  third-order  accurate  piecewise  spline  interpolation 
of  the  function  £(z).  As  one  can  see  in  Fig.4.10,  these  optimal  grid  results  exhibit  the  convergence  rate  of 
essentially  0(Af4)  and  provide  higher  accuracy  than  those  calculated  on  the  optimal  grid  Eq.(2.53),(2.54). 

To  demonstrate  the  superiority  of  the  optimal  grid  adaptation  over  the  equispaced  grid  point  distribution 
an  error  convergence  plot  of  a  symmetric  fourth-order  accurate  approximation  of  fx  calculated  on  a  uniform 
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Fig.  4.9.  Error  convergence  of  a  second-order  hybrid  approximation  calculated  with  the  consistent  and  inconsistent  dis¬ 
cretizations  of  the  metric  coefficient  on  the  optimal  and  uniform  grids. 


grid  with  the  same  number  of  grid  points  is  also  depicted  in  Fig.4.10.  The  L2  norm  of  the  truncation  error 
of  the  third-order  approximation  Eq.(4.10)  on  the  optimal  grid  is  reduced  by  a  factor  of  several  hundred  in 
comparison  with  the  fourth-order  accurate  results  obtained  on  the  uniform  grid. 

Error  convergence  plots  for  the  following  hybrid  approximation 

fdf\  f  gSf  (-2/i-i  -  3/i  +  6/f+i  -  fi+2) ,  i  even 
’  \8iJi  \  (-11  fi  + 18 fi+1- 9 fi+2 +  2 fi+3),  i  odd 

calculated  on  the  optimal  and  corresponding  uniform  grids  are  shown  in  Fig.4.11.  The  optimal  grid  results 
are  about  4-5  orders  of  magnitude  more  accurate  than  those  obtained  on  the  finest  uniform  grid.  However,  if 
the  metric  coefficient  is  evaluated  by  Eq.(4.10)  in  the  entire  computational  domain  while  the  approximation 
of  remains  the  same  Eq.(4.11)  the  error  convergence  rate  of  this  inconsistent  discretization  becomes  even 
less  than  0( A£3)  as  the  grid  is  refined. 

The  next  test  example  is  a  fourth-order  accurate  approximation  of  the  first  derivative  of  the  function 
f(x)  =  xm ,  where  the  parameter  m  has  been  chosen  to  be  49.  The  first  derivatives  and  x $  are  discretized 
by  a  five-point  symmetric  approximation 

(4.12)  (gt)i  =  Y2A£  ^ 9i~ 2  ~~  i_1  +  “  9i+2^ 5 

where  p(£)  is  either  /(£)  or  #(£).  It  can  be  shown  that  if  the  order  of  approximation  p  is  an  even  number 
then  for  /( x)  =  x™  the  optimal  grid  generation  equation  Eq.(2.14)  can  be  solved  analytically.  Thus,  we  have 

(4.13)  %opt{0  - 

The  above  mapping  is  optimal  in  the  sense  of  the  minimization  of  the  leading  truncation  error  term  if  m>  p 
otherwise  any  pth-order  accurate  difference  expression  approximates  the  first  derivative  fx  exactly.  If  we  fix 


27 


-0 -  Uniform,  3rd  order 


-m -  Numerical  optimal 


Fig.  4.10.  Error  convergence  for  a  third-order  approximation  of  fx ,  f(x)  =  xm  calculated  on:  1)  uniform  grid,  2)  optimal 
grid  generated  numerically,  3)  analytical  optimal  grid,  4)  grid  adapted  in  accordance  with  the  arc  length  criterion,  5)  uniform 
grid  using  fourth-order  accurate  discretization. 


the  parameter  m  to  be  sufficiently  large  one  can  observe  that  as  the  order  of  approximation  p  is  increased 
the  optimal  grid  Eq.(4.13)  becomes  more  uniform  that  correlates  with  the  above  theoretical  analysis.  The 
optimal  grid  point  distribution  can  also  be  calculated  numerically  by  using  Eq.(2.54).  At  each  grid  point 
the  unknown  parameter  a(x)  is  found  as  a  solution  of  the  equation 


(4.14) 


T4(a)  =  0, 


where  T4(a)  is  given  by  Eq.(2.59).  For  this  particular  choice  of  the  function  f(x ),  Eq.(4.14)  can  be  solved 
analytically  that  yields 


(4.15) 


1  777  -  4 
5  777  —  2  ' 


Note  that  the  optimal  mapping  Eq. (2. 54), (4. 15)  is  identical  to  Eq.(4.13)  if  we  set  p  =  4  in  it.  Error 
convergence  plots  calculated  on  the  analytical  Eq.(4.13)  and  numerical  Eq. (2. 54), (4. 15)  optimal  grids  as  well 
as  on  the  corresponding  uniform  grid  are  shown  in  Fig.4.12.  As  one  can  see  in  this  figure  the  fourth-order 
approximation  Eq.(4.12)  on  the  optimal  grid  Eq.(4.13)  exhibits  even  a  higher  convergence  rate  than  0(A£5) 
that  allows  one  to  reduce  the  L2  norm  of  the  truncation  error  by  6  orders  of  magnitude  compared  to  the 
uniform  grid  results.  The  numerical  approximation  of  both  the  second  derivative  and  the  integral  in  Eq.(2.54) 
leads  to  that  the  optimal  grid  Eq.(2.54),(4.15)  generated  numerically  provides  super  convergent  results  only 
on  coarse  grids  while  as  the  grid  is  refined  the  order  of  approximation  deteriorates  to  4.  Nevertheless,  the 
evaluation  of  fx  on  the  80-cell  optimal  grid  Eq.(2.54),(4.15)  is  about  3  orders  of  magnitude  more  accurate 
than  that  on  the  uniform  grid  with  the  same  number  of  grid  points.  One  of  the  main  reasons  of  such  a 
behavior  is  an  error  introduced  by  the  numerical  approximation  of  fxx  in  Eq.(2.54).  As  mentioned  above, 
the  optimal  mapping  Eq.(4.13)  is  singular  at  £  —  0  that  considerably  decreases  the  accuracy  of  the  fxx 
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Fig.  4.11.  Error  convergence  of  a  third-order  hybrid  approximation  calculated  with  the  consistent  and  inconsistent  dis¬ 
cretizations  of  the  metric  coefficient  on  the  optimal  and  uniform  grids. 


approximation  Eq.(4.1).  This  perturbation  introduced  into  the  optimal  grid  by  the  numerical  evaluation 
Eq.(4.1)  destroys  the  superconvergence  property.  However,  if  one  uses  the  exact  expression  for  fxx  despite 
that  the  integral  in  Eq.(2.54)  and  x(£)  are  calculated  numerically,  the  order  of  approximation  is  practically 
recovered  to  its  optimal  value  that  can  be  seen  in  Fig. 4. 12. 

To  demonstrate  the  importance  of  the  consistent  approximation  of  and  error  convergence  plots 
calculated  using  different  hybrid  approximations  on  the  optimal  and  corresponding  uniform  grids  are  depicted 
in  Fig. 4. 13.  The  fourth-order  accurate  hybrid  approximation  is  constructed  as  follows 


(4.16) 


Yskf  (fi-2  ”  8/i-l  +  8/i+l  -  /»+ 2)  > 

12^  (— 3/i-i  -  10/i  +  18/i+i  -  6 /i+2  +  A+s)  , 


i  even 
i  odd 


If  the  metric  coefficient  is  evaluated  by  the  same  difference  expression  employed  for  the  first  derivative 
Eq.(4.16)  then  the  leading  term  of  the  truncation  error  is  vanished  on  the  optimal  grid  Eq.(4.13).  It  is 
evident  in  Fig.4.13  that  the  truncation  error  of  the  consistent  hybrid  approximation  of  /$  and  exhibits 
the  convergence  rate  of  0(A£5).  At  the  same  time,  if  the  metric  coefficient  is  discretized  by  the  symmetric 
fourth-order  accurate  formula  Eq.(4.12)  in  the  entire  computational  domain,  while  the  same  approximation 
Eq.(4.16)  is  used  for  /€,  the  convergence  rate  deteriorates  to  0(A£4)  and  the  truncation  error  increases  by 
a  factor  of  50-100  in  comparison  with  the  consistent  approximation  results.  The  deterioration  of  the  error 
convergence  rate  on  the  finest  optimal  mesh  is  presumably  caused  by  the  machine  accuracy. 


4.2.  2D  Test  Example.  We  shall  seek  a  particular  solution  of  Eq.(3.6)  in  the  following  form 


(4.17) 


/&»?)= 

Xopti^rj)  = 

yoPt(^r])  =  e0^\ 
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Fig.  4.12.  Error  convergence  for  a  fourth-order  approximation  of  fx>  /(#)  =  xm  calculated  on:  1)  uniform  grid,  2) 
optimal  grid  generated  numerically,  3)  analytical  optimal  grid,  4)  grid  adapted  in  accordance  with  the  arc  length  criterion,  5) 
uniform  grid  using  fifth-order  accurate  discretization, 6)  numerical  optimal  grid  generated  with  the  exact  fXx- 


where  a ,  /?,  and  7,  </>,  9,  ip  are  given  and  unknown  constants,  respectively.  Note  that  this  choice  of  /, 
x ,  and  y  uniquely  defines  the  function  f(x,y)  in  the  physical  domain.  Since  the  above  mapping  must  be 
nonsingular  the  Jacobian  of  the  mapping,  which  is 

(4.18)  =  (rl>  -  ^)e(7+*)?e<*+^, 
should  be  positive  in  the  whole  computational  domain  that  leads  to 

(4.19)  7 ip-(p9>  0. 

Substituting  Eq.(4.17)  into  the  first  equation  of  Eq.(3.6)  yields 

(4.20)  (7^  “  <p0)a3  =  (—(pa  +  7 0)93  +  (ipa  -  9(5) 7s. 

Equation  (4.20)  together  with  the  constraint  Eq.(4.19)  give  us  a  family  of  the  optimal  grids.  The  equation 
is  simplified  considerably  if  we  assume  that  <p  —  ip  =  f3  —  1.  Under  this  assumption  Eq.(4.20)  and  (4.19)  are 
reduced  to 

(4.21)  (7  -  #)a3  =  (7  -  oi)93  4-  (a  -  9) 7s 

(4.22)  7  “  0  >  0, 
correspondingly.  Equation  (4.21)  has  three  real  roots 

7i  —  a  —  9 

(4.23)  72  =  0 

73  =  a 
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Fig.  4.13.  Error  convergence  of  a  fourth-order  hybrid  approximation  calculated  with  the  consistent  and  inconsistent 
discretizations  of  the  metric  coefficient  on  the  optimal  and  uniform  grids. 

The  roots  72  and  73  are  not  appropriate  because  the  second  root  does  not  meet  the  inequality  Eq.(4.22) 
while  the  third  root  implies  that  f(x)  =  x.  Therefore,  the  only  non-trivial  solution  of  Eq.(4.21),  (4.22)  is 
7  +  6  =  a.  Introducing  a  parameter  m  so  that  7 /Q  —  m  the  particular  solution  of  Eq.(3.6)  can  be  written  in 
the  following  form 

v  -ma  £ 

Xopt(£,r))  =e"'+i*e7' 

(4.24)  yopt(£,v)  = 

m+2  2m  +1 

f(X,y )  =  X  rn-lym-l  . 

In  the  present  test  example  the  parameters  m  and  a  have  been  chosen  to  be  10  and  3,  respectively.  The 
corresponding  optimal  41  x  21  grid  and  30  isolines  of  the  function  f(x,y)  are  depicted  in  Fig. 4. 14.  It  is 
notable  that  the  optimal  grid  is  orthogonal  neither  in  the  domain  nor  at  the  boundaries.  Moreover,  the  grid 
lines  are  concentrated  near  strong  gradients  and  at  the  same  time,  they  are  not  strictly  aligned  to  the  isolines 
of  f(x,  y).  A  second-order  accurate  approximation  of  fx  is  obtained  by  using  two-point  central  differences  for 
all  the  derivatives  in  Eq.(3.1).  A  uniform  grid  is  generated  by  the  transfinite  interpolation  of  the  boundary 
nodes,  which  are  uniformly  distributed  along  the  boundaries.  Since  the  optimal  grid  Eq.(4.24)  has  been 
constructed  under  the  assumption  that  the  leading  term  of  the  truncation  error  in  the  f  coordinate  vanishes 
on  the  optimal  grid  we  shall  refine  the  grid  only  in  £  while  the  number  of  grid  cells  in  7]  is  fixed  and  equal 
to  20.  Note  that  the  grid  refinement  in  the  7]  coordinate  makes  no  influence  on  the  convergence  rate  of  the 
truncation  error  that  is  consistent  with  Eq.(3.6). 

A  comparison  of  the  truncation  error  convergence  obtained  on  the  optimal  and  uniform  grids  is  shown  in 
Fig.4.15.  Similarly  to  the  ID  test  examples,  the  global  order  of  the  symmetric  second-order  approximation 
in  two  dimensions  is  increased  by  more  than  one  on  the  optimal  grid.  Furthermore,  the  norm  of  the 
truncation  error  is  about  4  orders  of  magnitude  less  than  that  obtained  on  the  corresponding  uniform  grid. 
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Fig.  4.14.  Optimal  40  x  20  grid  and  30  isolines  of  the  function  f(x). 


As  can  be  seen  in  Fig.4.15,  the  new  grid  adaptation  criterion  enables  one  to  reach  the  asymptotic  convergence 
rate  on  coarse  grids  while  the  application  of  a  third-order  accurate  discretization  on  the  uniform  grid  does 
not  permit  us  to  get  so  essential  reduction  in  the  truncation  error  as  on  the  optimal  grid. 

The  importance  of  the  identical  approximation  of  the  first  derivatives  and  fv  and  the  metric  coefficients 
y and  xn,  yn ,  respectively  is  illustrated  in  Fig.4.16.  The  figure  shows  that  if  /$,  ,  and  y $  are  evaluated 

by  the  same  hybrid  discretization  Eq.(4.3)  the  order  of  approximation  in  £  is  increased  by  one  if  grid  points 
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Fig.  4.15.  Error  convergence  for  a  second-order  approximation  of  fx  calculated  on:  1)  uniform  grid ,  2)  analytical  optimal 
grid,  and  3)  uniform  grid  using  third-order  accurate  discretization. 


are  redistributed  in  accordance  with  Eq.(4.24)  regardless  what  second-order  approximations  are  used  for 
Xq,  and  Ur).  However,  if  the  metric  coefficients  x f  and  are  evaluated  by  the  two-point  symmetric  second- 
order  difference  expression  in  the  entire  computational  domain,  whereas  both  the  hybrid  approximation  of  f f 
Eq.(4.3)  and  the  optimal  grid  Eq.(4.24)  remain  the  same,  the  order  of  approximation  of  fx  in  f  deteriorates 
to  2  and  the  truncation  error  is  increased  by  a  factor  of  103. 

5.  Conclusion.  The  new  grid  adaptation  strategy  based  on  the  minimization  of  the  leading  truncation 
error  term  of  an  arbitrary  pth-order  finite  difference  discretization  has  been  developed.  The  basic  idea  of  the 
method  is  to  redistribute  grid  points  so  that  the  leading  truncation  error  terms  due  to  the  differential  operator 
and  the  metric  coefficients  cancel  each  other  so  that  the  design  order  of  approximation  on  the  optimal  grid 
is  increased  by  one  in  the  entire  computational  domain.  In  contrast  to  most  of  the  adaptive  grid  techniques, 
for  the  present  method  neither  the  truncation  error  estimate  nor  the  specification  of  weighting  parameters  is 
required.  Another  very  attractive  characteristic  of  the  new  approach  is  its  applicability  to  hybrid  discretiza¬ 
tions.  It  has  been  proved  that  if  the  differential  operator  and  the  metric  coefficients  are  evaluated  identically 
then  the  same  optimal  grid  adaptation  criterion,  which  is  valid  for  non-hybrid  discretizations,  can  be  used 
in  the  entire  computational  domain  regardless  of  points  where  the  hybrid  difference  operator  switches  from 
one  approximation  to  another.  One  of  the  main  advantages  of  the  new  method  is  that  it  can  be  directly 
extended  to  multiple  dimensions.  It  has  been  shown  that  the  multidimensional  grid  adaptation  criteria  are 
fully  consistent  with  the  one-dimensional  counterpart.  The  ID  and  2D  numerical  calculations  show  that 
the  truncation  error  obtained  on  the  optimal  grid  is  both  superconvergent  and  reduced  by  several  orders 
of  magnitude  in  comparison  with  the  uniform  and  standard  adaptive  grid  results  for  all  the  test  examples 
considered. 
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Fig.  4.16.  Error  convergence  of  a  second-order  hybrid  approximation  calculated  with  the  consistent  and  inconsistent 
discretizations  of  the  metric  coefficient  on  the  optimal  and  uniform  grids. 
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